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ABSTRACT 

A  formulation  of  the  equations  governing  the  large  transient  defor¬ 
mation  of  Kirchhoff  shells  especially  suitable  for  use  in  finite  differ¬ 
ence  structural  response  codes  is  presented.  It  is  shown  how  these 
equations  are  employed  in  REPSIL,  a  computer  code  developed  at  foe  BRL 
for  predicting  foe  large  deflection  of  thin  elastoplastic  shells  subjected 
to  blast  and  iwpuise  loadings .  The  REPSIL  formulation  is  cohered  with 
foe  somewhat  similar  formulations  used  by  PETROS  1  and  PETROS  2,  two 
closely  related  codes  developed  by  MIT. 
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I .  INTRODUCTION 

This  report  is  mainly  concerned  with  presenting  c  formulation  of 
the  equations  governing  'J:c  large  transient  deformations  of  shells,  es¬ 
pecially  suitable  ror  use  in  finite  difference  computer  codes.  The 
recent  adaptation  of  these  equations  to  REPSIL,  a  code  developed  at  the 
BRL  for  predicting  the  large  deflections  of  thin  shells  subjected  to 
blast  and  impulse  loadings,  has  resulted  in  a  number  of  improvements, 
which  are  also  described  in  this  report.  While  the  equations  as  formu¬ 
lated  in  the  main  body  of  the  report  are  quite  general,  being  restricted 
only  by  the  Kirchhoff  hypothesis,  as  applied  in  REPSIL  they  are  further 
restricted  to  thin  elastoplastic  shells  for  which  the  rotatory  inertia 
is  negligible.  A  detailed  description  of  the  REPSIL  code  will  be  in¬ 
cluded  in  a  forthcoming  user's  manual. 


1.1  Historic  Background 

In  order  to  put  this  report  and  the  associated  computer  code  REPSIL 
in  their  proper  relation  to  work  done  at  the  Aeroelastic  and  Structures 
Research  Laboratory  of  MIT  on  similar  codes,  a  short  history  cf  the  de¬ 
velopment  of  REPSIL  is  in  order.  REPSIL  has  evolved  from  a  computer 

♦ 

code  developed  by  LEECH  [1]  originally  called  PETROS  and  now  called 
PETROS  1  since  the  recent  development  by  MORINO,  LEECH  and  WITMER  [2]  of 
the  improved  PETROS  2.  For  its  time  PETROS  1  was  the  most  sophisticated 
code  available  treating  the  large  transient  response  of  shells.  For  this 
reason  it  was  one  of  the  codes  selected  for  use  in  the  numerical  analysis 
portion  of  the  research  program  on  shells  and  structures  at  the  BRL. 

Upon  close  study,  BUFFINGTON  [3]  found  that  PETROS  1  contained 
a  number  of  flaws,  in  addition  to  having  a  rather  limited  choice  of 
printing  options.  He  discovered  that  the  clamped  edge  boundary  condition 


Numbers  in  brackets  refer  to  the  list  of  references  on  page  123. 


was  incorrectly  formulated  and  that  the  linear  approximation  employed  in 
the  stress  increment  calculation  during  plastic  flow  gave  stress  states 
outside  the  yield  surface.  He  revised  the  code,  by  increasing  the  print¬ 
ing  options,  adding  plotting  capability,  correcting  the  clamped  edge  con¬ 
dition  and  reformulating  the  stress  calculation  based  on  the  exact  quad¬ 
ratic  expression.  The  revised  code  was  named  REPSIL.  The  details  of 

* 

the  revised  stress  calculation  are  documented  in  [4]  . 

In  this  early  version,  however,  REPSIL  still  retained  many  of  the 
uisatisfactory  features  of  PETROS  1.  For  example,  only  cylindrical 
panels  with  clamped  edges  could  be  treated,  the  finite  difference  mesh 
had  to  be  square,  the  strain-displacement  relation  did  not  properly 
account  for  the  curvature  cf  the  shell,  the  finite  difference  derivatives 
along  the  boundaries  were  inconsistent,  etc.  Consequently,  a  complete 
reformulation  of  REPSIL,  including  its  theoretical  foundation,  was  begun 
at  the  BRL.  The  reformulation  has  advanced  sufficiently  to  warrant 
documentation  of  the  present  version  of  REPSIL.  This  report  documents 
the  mathematical  model  on  which  REPSIL  is  presently  based.  A  user's 
manual  documenting  the  computational  reformulation  of  REPSIL  will  follow 
shortly. 

The  Aeroelastic  and  Structures  Research  Laboratory  of  MIT  has  also 
undertaken  a  complete  revision  of  PETROS  1,  which  has  resulted  in  the 
code  PETROS  2  [2],  Although  REPSIL  and  PETROS  2  share  a  number  of 
features,  they  are  independent  extentions  of  PETROS  1  and  do  not  coincide. 
The  differences  between  REPSIL  and  PETROS  2,  as  well  as  their  differences 
from  PETROS  1,  are  discussed  in  Chapter  8. 


I 


PETROS  2  [2]  uses  a  slightly  modified  version  of  this  stress  calculation. 

tile 

My  colleague.  Dr.  Huffington,  was  responsible  for  first  making  me 
aware  of  some  of  these  flaws. 
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1.2  Outline  of  Remainder  of  Report 


Tne  next  five  chapters.  Chapters  2  through  6,  develop  rigorously 
the  theory  of  the  motion  of  shells  subject  to  the  Kirchhoif  hypothesis. 
Although  a  functional  relation  is  presumed  to  exist  between  the  stress 
and  strain  components,  none  is  explicitly  used  to  develop  the  theory, 
in  keeping  with  the  general  nature  of  these  chapters. 

Chapter  7  is  more  directly  concerned  with  REPSIL.  First  the  general 
theory  developed  previously  is  specialized  to  thin  shells  and  the  rota¬ 
tory  inertia  is  neglected.  Moreover,  the  stress  -  strain  relation  for 
any  material  conposing  the  shell  is  assumed  to  be  isotropic  and  linearly 
elastic  -  perfectly  plastic.  These  assumptions  give  us  a  convenient 
version  of  the  equations  forming  the  mathematical  basis  of  REPSIL,  which 
are  then  used  to  illustrate  the  computational  algorithm  employed  by 
REPSIL  to  solve  problems.* 

As  mentioned  earlier.  Chapter  8  gives  a  comparison  of  REPSIL, 

PETROS  1  and  PETROS  2.  First  the  theoretical  and  then  the  computational 
differences  are  discussed.  These  are  then  summarized  in  tabular  form. 

The  reader  principally  interested  in  Chapters  7  and  8  is  advised 
to  read  at  the  very  minimum  Chapter  2  in  order  to  acquaint  himself  with 
the  differential  geometric  approach  used  in  REPSIL,  as  well  as  in 
PETROS  1  and  PETROS  2,  and  if  further  time  permits,  he  should  skim 
Chapters  3  through  5  to  familiarize  himself  with  the  significance  of 
the  more  commonly  encountered  terms.  Chapter  6  can  be  safely  skipped. 


Although  the  algorithm  is  illustrated  assuming  perfectly  plastic  ma¬ 
terial  behavior,  actually  REPSIL  has  the  capability  of  treating  strain 
hardening,  strain  rate  sensitive  materials  through  the  use  of  the 
mechanical  sublayer  model  described  in  [2;  Sect.  5.4.2]. 
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2.  GEOMETRY  OF  A  DEFORMING  SHELL 

In  this  chapter  we  present  the  mathematics  needed  to  describe  the 
motion  of  a  shell  geometrically.  This  involves  introducing  cuncepts 
and  results  from  differential  geometry.  We  first  treat  the  motion  of 
the  reference  surface  of  the  shell,  after  which  the  consequences  of 
assuming  the  Kirchhoff  hypothesis  for  the  deformation  through  the  shell 
thickness  are  obtained. 


2.1  Motion  of  the  Reference  Surface 


A  shell  can  be  described  qualitatively  as  a  three-dimensioned  body 

having  one  of  its  dimensions,  the  thickness,  small  in  comparison  to  its 

other  dimensions.  Shell  theory  takes  advantage  of  this  property  in 

order  to  simplify  the  studty  of  shells  by  reducing  an  essentially  three- 

dimensional  theory  to  one  defined  on  a  two-dimensional  subspace,  i.e.  a 

surface.  This  is  accomplished  by  assuming  that  the  deformation  through 

the  thickness  is  of  a  known  functional  form.  Once  this  is  done,  we  need 

only  study  the  deformation  of  a  given  surface  enfcedded  in  the  shell, 

* 

called  the  reference  surface  ,and  the  associated  changes  in  the  shell 
variables  defined  on  this  surface. 

a** 

Denoting  the  surface  parameters  by  £  ,  the  time  by  t  and  the 

position  vector  in  Euclidean  3-space  relative  to  some  fixed  origin  by 


i. 


I 


It  is  customary  to  pick  the  reference  surface  to  be  the  middle  surface 
of  the  shell;  that  is,  the  surface  equidistant  From  the  bounding  faces 
of  the  shell.  However,  whether  a  surface  initially  equidistant  from 
the  bounding  faces  remains  a  middle  surface  throughout  the  history  of 
the  deformation  is  a  function  of  the  form  assumed  for  the  deformation 
through  tne  thickness.  Hence,  when  no  form  has  been  explicitly  assumed, 
the  term  reference  surface  is  to  be  preferred. 

Index  notation  is  employed  throughout  this  report,  with  Greek  indices 
being  tacitly  assumed  to  range  over  the  integers  1,  2.  Hence, 

5“  s  {>,  C2  • 
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r  ,  the  motion  of  the  reference  surface  is  specified  by  the  reference 
surface  deformation  function: 


r  =  r  U  ,t)  , 


(2.1) 


with  £  ranging  over  some  connected  region  in  the  plane  and  t  over 
some  interval  of  the  real  line.  We  demand  that  the  mapping  £a  -*■  r  be 
one-to-one  for  any  time  t.  We  also  require  that  the  function  r  (£a,t) 
be  sufficiently  differentiable  so  as  to  guarantee  the  continuity  of  all 
derivatives  subsequently  introduced.  Moreover,  we  assume  for  simplicity 
that  the  £a  are  material  coordinates;  that  is,  for  fixed  fa  the  image  r 
of  the  point  £a  under  the  mapping  (2.1)  represents  the  successive  po- 
siions  of  the  same  material  point  or  particle  on  the  reference  surface. 


2.2  Differential  Geometry  of  the  Reference  Surface 


In  this  section  we  are  concerned  with  the  spatial  characterization 
of  the  reference  surface  at  some  fixed  instant  of  time.  We  will  need 


tc  use  the  language  of  the  differential  geometry  of  surfaces  in  Euclidean 
3-space.  Our  treatment  of  this  subject  is  not  meant  to  be  exhaustive.* 
Rather,  we  aim  to  present  only  selected  results  with  a  minimi*  of  proof. 
In  this  section  we  shall  not  bother  to  indicate  the  dependence  of  vari¬ 
ables  on  t  explicitly. 

By  holding  the  parameter  S1  constant  and  allowing  to  vary,  we 
generate  a  member  of  the  family  of  =  const,  coordinate  curves. 
Conversely,  by  holding  S2  constant  and  varying  we  generate  a  member 
of  the  family  of  S2  =  const,  coordinate  curves.  We  require  that  these 
two  families  of  curves  be  independent  in  the  sense  that  no  member  of 
one  be  anywhere  tangent  to  any  member  of  the  other.  Thus,  these  two 
families  define  a  curvilinear  coordinate  system  on  the  reference  surface. 


Extended  treatments  of  the  differential  geometry  of  surfaces  can  be 
fouid  in  numerous  text  books  on  tensor  analysis  and  differential 
geometry,  of  which  the  following  are  a  representative  collection: 
EISENHART  [6],  SOKOLNIKOFF  [7]  and  THOMAS  [8]. 
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The  partial  derivative  of  the  function  r  (5°)  with  respect  to  C1 
holding  €2  constant 

h  r  >x  (cV  (2-2) 

defines  a  vector  tangent  to  a  £2  =  const  coordinate  curve  and  the  other 
partial  derivative  holding  £,x  constant 

a2  U“)  =  r  ,2  (C°)  (2.3) 

a  vector  tangent  to  a  S1  =  const  coordinate  curve.  We  require  that 
the  vectors  a  nowhere  vanish.  Hence,  by  virtue  of  the  previously 
assumed  independence  of  the  families  of  coordinate  curves,  the  pair  aQ 
form  a  linearly  independent  set  of  vectors  tangent  to  the  reference 
surface  and  hence  serve  as  a  basis  for  vectors  tangent  to  the  reference 
surface.  They  are  in  fact  the  basis  generated  by  the  coordinate  curves. 

The  unit  normal  to  the  surface  is  sinply 

„  ai  x  ao 

n  =  A - .  (2.4) 

Uj  ^  a2| 

Clearly  the  triad  (ao>  n)  forms  a  basis  for  vectors  in  3-space,  with 
the  special  orthogonality  property  that 

n  •  a  =  0  .  (2.5) 

-a 


* 

We  use  the  comma  notation  for  partial  derivatives  with  respect  to 
material  coordinates:  a  partial  derivative  with  respect  to  £a  is 
denoted  by  a  subscripted  comma  followed  by  the  index  a  ;  that  is. 


3 
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The  covaricnt  components  of  the  metric  or  the  first  fundamental 
terror  of  the  reference  surface  are  defined  as 


a  6 


=  a 
_e 


(2.6) 


They  form  a  symmetric  and  positive  matrix.  Denoting  the  determinant  of 
the  metric  matrix  by  a,  we  have 


a  =  Det  aag  =  a^  -  (a12)2  =  x  a2|2  >  0  .  (2.7) 

The  con iravari an t  components  of  the  surface  metric  can  then  be  defined 
as 


a11  =  »22^a  »  &12  =  ~  al2//a  =  ^  ‘  &22  =  aH^a  •  (2-8) 

and  they  also  form  a  symmetric  and  positive  matrix.  Moreover,  a  P  is 
the  inverse  of  a 


(2-9) 


with  (=  1  when  a  =  6,  =0  when  a  /  B)  the  Kronecker  delta.  By  era- 
B 

ploying  the  contravariant  components  of  the  metric  tensor,  the  dual 
basis  is  generated: 


a  aB 
a  =  a  ag  . 


(2.10) 


The  dual  is  inversely  related  to  the  basis  a^  and  is  also  normal  to  the 


We  employ  the  summation  convention:  terms  having  a  repeated  index,  once 
as  a  subscript  and  once  as  a  superscript,  are  to  be  summed  over  the 
range  of  that  index. 
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wiiiyjLw#^^ 


The  covariant  and  contravariant  components  of  the  skew-sumnetr-ic 
* 

surface  tensor  z  are  defined  as 


eaB  “  &2  ®aB 


(2.12) 


where  e  _  =  e  (=1  when  0=153=2,=  -1  when  a  =  2  5  S  =  1, 
afc5 

=  0  when  a  =  B)  is  the  permutation  symbol.  The  components  of  e  and  the 
Kronecker  delta  satisfy  the  usual  relations: 


Y  5 


5  Y 


(2.13) 


Using  the  e  tensor,  we  can  summarize  the  relationships  between  the  basis 
a.Q  ,  the  dual  basis  aa  and  the  surface  normal  n  in  the  following  compact 
forms: 

v  a  B  ag  ,,  ... 

a  *  a„  =  e  „n  ,  a  x  a  =  e  n  ,  (2.14) 

_a  _B  ap_  _  _ 

B  a  aB 

nxa=eQa  ,  n  x  a  =  e  aQ. 

_a  aB-  *  _B 

A  tensor  as  important  as  the  metric  or  first  fundamental  tensor  for 
surfaces  embedded  in  3-space  is  the  second  fundamental  tensor ,  which 
measures  the.  spatial  rate  of  change  of  the  normal  n  with  respect  to  the 
surface  basis.  Its  covariant  components  are  defined  as 


b«e  ■  -  ?« •  ;-e =  ;  •  • 


(2-15) 


5  is  only  a  tensor  under  coordinate  transformations  with  positive 
Jacobians  since  the  permutation  symbol  eag  transforms  like  a  relative 
tensor  of  weight  ±1;  cf.  SOKCLN'IKOFF  [7;  p.  107]. 
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the  last  equality  clearly  showing  that  the  b  tensor  is  symmetric.  Simi¬ 
larly,  the  components  of  the  connection  for  the  surface  are  defined  as 


(2.16) 


with  the  symmetry  with  respect  to  the  first  and  second  indices  exhibited 

★ 

by  the  last  equality.  A  standard  reduction  relates  the  components  of 
the  connection  to  the  partials  of  the  metric  components 


1  y6 

2 


^aB6,a  +  a6a,B 


(2.17) 


Hence,  the  components  of  the  connection  are  the  Christoffel  symbols. 
It  easily  follows  from  this  that 


% 


a  „  =  T  a 
,a  2 


1  Jz  By 


o-o  =  a 2  r  ^ 

BYjCi  aB 


(2.18) 


The  covariant  derivative  of  any  quantity  with  respect  to  £a,  denoted  by 
an  attached  subscript  I o,  is  defined  in  the  usual  way  through  the 
connection.  The  surface  metric  and  the  surface  tensor  e  have  the 
property  that  the  covariant  derivative  of  their  components  identically 
vanish: 


aaB  lY 


=  0 


*  GaB  I Y 


=  0 


eaBjY  =  0.  (2.19) 


See  SOKOLNIKOFF  [7;  p.  129]. 


**  3 

For  example,  the  covariant  derivative  of  the  two  index  quantity  £ 
with  respect  to  Ca  is 


& 


B 


M6  r  6 


Y  I  o  °Y,a  0  y  d5‘oy  * 

cf.  EISENHART  [6;  p.  110]  or  THOMAS  [8;  p.54]. 
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r  ~  .-*53^.?^.. 


This  well  kncvm  property  shall  be  freely  used  throughout.  It  easily 
follows  from  this  property  and  the  definition  of  the  second  fundamental 
tensor  that 


?a  I  6  =  baB  *  ’  ?°l  B  “  b6  -  *  "  I  6  ~  "  bS  *a 


(2.20) 


Note  that  the  components  of  the  metric  are  used  to  raise  and  lower  the 
indices  of  the  second  fundamental  tensor  (i.e.  =  a°Y  b^),  as  is 

commonly  done  with  the  components  of  any  tensor  associated  with  the 
surface.  Lastly,  we  make  mention  of  the  fact  that  due  to  the  second 
fundamental  tensor  and  the  connection  being  defined  as  derivatives  of 

7t 

the  basis  a^  they  satisfy  the  Codazzi  and  Gauss  equations 

b  -  b  .  .  =  0  ,  R  o  *  =  b  .b.  -  b  bot  ,  (2.21) 

aBI  y  ayl  B  aBy6  a<5  By  ay  B6  v 


where 


R«M  ■  V,  -  V'Y  *  V  V  -  V  rvvV)  (2'22) 


is  the  Riemann  Chris toffel  curvature  tensor  for  the  surface. 


2.3  Kinematics  of  the  Reference  Surface 

The  time  derivative  holding  the  material  coordinates  constant  is 
the  material  derivative .  The  material  derivative  of  the  function 
r  (5a,t)  is  the  velocity  at  the  time  t  of  the  particle  with  material 
coordinates  5°: 


v  (e  ,t)  =  r  U“,t)  . 


(2.23) 


Resolved  into  normal  and  tangential  components,  the  velocity  becomes 


v  =  v  a  +  v  n  - 
_a 


(2.24) 


Cf.  EISENHART  [6;  p.  219]  or  THOMAS  [8;  p.  89]. 

k 

A  superposed  dot  denotes  the  material  derivative. 


The  covariant  derivative  of  the  velocity,  the  velocity  gradient 3  is  also 
resolved  into  normal  and  tangential  components 


v,  =  v  =  (v6.  -  b^  v)  a  +  (v.  +  b  0v^)  n  , 
,la  ,,a  v  la  a  1  .6  v  la  aS  '  , 


(2.2S) 


using  relations  (2.20).  A  reversal  of  order  of  differentiation  shows 
that  the  velocity  gradient  is  the  material  derivative  of  the  basis: 


v,  =  k 
,  I  a  .a 


(2.26) 


The  divergence  of  the  velocity  v  can  be  written  as 


„  _  a  a  ,  a 

V-v  =  a  *  v.  =  v  .  -b  v  . 
,  -  .  I  a  la  a 


(2.27) 


Combining  the  last  two  equations,  we  obtain  the  useful  relation 


'h  1  k  aB  *  h  „ 

a  =  2  3  3  aaB  =  3  V*! 


(2.28) 


Since  the  normal  vector  n  is  of  constant  (unit)  magnitude,  its 
material  derivative  w  will  be  tangential  to  the  surface  and  hence  have 
the  resolution: 


.  a  a 

n  =  cj  =  (o  a  =  u>  a 
-  .a  a  _ 


(2.29) 


Also,  since  n  and  aa  are  normal,  their  material  derivatives  cannot  be 
totally  independent,  but  must,  indeed,  satisfy  the  relation 


u*a  +  n  •  v  ,  =  0 

,a  _  ,  la  » 


(2.3C) 


from  which  it  follows  that 


(o  =  -  n  •  v  .  a 
>  I  a  _ 


(2.31) 


or,  using  (2.25)  and  (2-29), 


w  =  -  (v  .  +  b  ,  v  )  . 

a  v  I  a  aB 


(2.32) 
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£ 


The  latter  allows  us  to  write  the  velocity  gradient  in  the  form 


8  8 

v.  =  (v  ,  -  bp  v)  a„  -  w  n 

Ja  '■ia  a  .8  ct  _ 


(2 . 33) 


Interesting  results  could  now  be  obtained  relating  the  material  deriva¬ 
tives  of  the  metric,  second  fundamental  tensor,  connection,  etc-  to  vari¬ 
ous  covaiiant  derivatives  of  the  velocity  components,  but  since  the  ma¬ 
terial  presented  suffices  for  subsequent  sections,  we  will  not  pursue 
this  line. 

2.4  Kirchhoff  Hypothesis 

The  Kirchhoff  hypothesis  assumes  that  the  shell  deforms  in  such  a 
way  that  particles  initially  on  a  normal  to  the  reference  surface  will 
find  themselves  after  the  shell  is  deformed  on  the  same  normal  with  their 
distance  from  the  reference  surface  unchanged.  This  hypothesis  is  an 
idealization  of  tne  observation  that  in  many  situations  very  little 
shearing  or  extension  occurs  through  the  thickness  over  major  portions 
of  the  shell  during  deformation. 


We  shall  assume  that  the  shell  is  of  uniform  thickness  and  identify 
the  reference  surface  with  the  middle  surface  of  the  shell,  although  any 
surface  parallel  to  the  middle  surface  would  give  an  equivalent  theory. 
Using  the  coordinate  system  introduced  previously,  we  have  by  the 
Kirchhoff  hypothesis  that  a  particle  located  a  distance  c,  from  the 
middle  surface  along  the  normal  n  passing  through  the  middle  surface 
paicicle  with  coordinates  £a  ,  will  always  be  situated  on  the  normal 
through  the  same  middle  surface  particle  at  the  same  distance  c.  Hence 
for  all  time  the  position  of  a  particle  is  specified  by  the  equation. 


x  =  r  (CC,t)  +  c  n  (Ca,t)  . 


(2.34) 


£ 

•2 

c 


i 


I 

<1 


.5 


as  indicated  in  Figure  2.1. 

* 

|  Notice  that  a  particle  off  the  middle  surface  will  always  have 

I  associated  with  it  the  same  coordinates  (5a,0,  with  c  J-  0.  Moreover, 

a  particle  on  the  middle  surface  is  always  associated  with  the  coordinates 

1 


3 


(Ca,0).  Hence  the  triplex  (C°,?)  defines  a  material  coordinate  system 
for  all  the  particles  composing  the  shell.  Taking  h  to  be  th  uniform 
thickness  dimension  of  the  shell,  s  will  range  between  the  limits  ±  y  , 
while  as  before  £a  range  over  some  suitable  connected  region  in  the  plane. 

Definitions  (2.23)  and  (2.29)  allow  us  to  write  the  particle  velocity 
as 

x  =  v  +  ?  to  (2.35) 

The  derivatives  of  (2.34)  with  respect  to  material  coordinates  (5°,?) 
define  a  basis  in  3-space: 

ix 

§B  ~  *  §3  =  3^  '  (2.36) 

From  (2.20)^  it  follows  that  the  basis  g^  is  related  to  the  basis 
(aQ,n)  through  the  relations 

?6  =  WB  -Y  *  §3  =  «  »  (2.37) 

where  the  convenient  abbreviation 

"2  =  5e  -  5  be  <2-38> 

is  employed.  The  pair  g„  (f  ,5)  defines  a  basis  for  the  lanella 
parallel  to  the  middle  surface  at  a  distance  £  .  When  5  =  0  the  basis 
gg  reduces  to  the  middle  surface  basis  ag  .  Notice  that  n  is  normal 
to  each  lamella,  since 

n  •  gQ  =0  .  (2.3S 


* 

Latin  indices  are  introduced  for  the  range  1,  2,  3. 
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A  metric  for  that  portion  of  space  occupied  by  the  shell  can  be 
specified  using  the  basis  as  follows 


gij 


=  §i 


(2.40) 


The  determinant  and  inverse  of  the  metric  are  then  defined  in  the  usual 
way 


1  ijk  lmn 

8  '  Det  *ij  '  6  e  6  8il  8jm  8kn 


ii  1  ike  jmn 

5  =  2g  6  eJ 


gkm  gln  * 


(2.41) 


x  i  K 

where  e  J  (=1  when  i,  j,  k  is  an  even  permutation  of  1,  2,  3;  =  -1 
when  i,  j,  k  is  an  odd  permutation;  =  0  when  i,  j,  k  not  distinct)  is 
the  3~dimensional  permutation  symbol.  The  inverse  allows  us  tc  define 
the  dual  basis 

g1  =  gij  gj  ,  (2.42) 

having  the  reciprocal  property 


i 

§ 


(2.43) 


The  metric  and  its  inverse  assume  the  simple  forms 


gij 


/hi  si2 

I  g21  g22 

\0  0 


(2.44) 


due  to  g^  being  a  unit  vector  normal  to  gQ  .  Hence  v:e  need  only  work 

with  the  metric  g  „  of  lamellar  parallel  to  the  middle  surface  and  its 

aBaB  i 

inverse  g  .  Moreover  the  determinant  g  and  dual  basis  g  become 
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n  1  ag  y6 

8  *  0et  Soe  ‘  2  e  e  8«y  gB6  - 


a  aB 

2  “  8  8g  » 


g3  =  g  =  n  . 
-3 


(2.45) 


Using  (2.37),  we  relate  the  above  quantities  to  the  corresponding 
middle  surface  quantities 


=  jr  « 


ot8  -la  -18  y6 

rr  -•  it  o 


ga8  =  VB  \6  »  *  =  w  8  U  6  3 


,  .2  a  -la  8 

g  =  (p)z  a  ,  g  =  u  B  a 


'’"Id  (X 

where  p  and  p  “  denote  the  determinant  and  inverse  of  p“  : 

p  8 


(2.46) 


1  air  8  6  -la  1  ay  6 

u=2e  eSi  ya  N  '  UB=P6  e86vv 


(2.47) 


Clearly,  the  basis  g  ,  its  dual  g  ,  the  metric  g  _  and  its 

_a  ^  ap 

inverse  only  exist  when  p  ^  0  .  It  can  be  shown  that  a  sufficient 
condition  to  guarantee  existence  for  c;  in  the  range  ±  y  is  h  <  2  R  .  , 

where  R.  is  the  minimum  principle  radius  of  curvature  of  the  middle 
surface.  Kith  this  condition  p  >  0  .  Ke  shall  assume  this  mild  con¬ 
dition  to  hold  in  developing  the  theory  of  general  Kirchhoff  shells,  in 
Chapters  2  through  6.  In  specializing  the  general  theory  to  that  forming 
the  basis  for  REPSIL  in  Chapter  7,  we  go  further  and  require  that  the 
more  stringent  thin  shell  hypothesis  holds,  wherein  h  «  R  .  .so  that 
terms  of  order  |s  b  |  can  be  neglected  in  conparison  with  unity. 

P 


See  NAGHDI  [9;  p  17]. 
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3.  STRAIN  IN  A  SHELL 

In  this  chapter  a  measure  of  strain  is  defined  based  on  the  change 
in  the  metric.  This  measure  is  simplified  using  the  Kirchhoff  hypothesis 
Equations  are  next  obtained  relating  the  incremental  change  in  strain 
to  the  gradients  of  the  displacement  increments  of  the  middle  surface. 


3. 1  Definition  of  Strain 

Strain  is  any  quantity  measuring  the  amount  of  stretching  that 
occurs  in  the  local  neighborhood  of  a  particle  during  deformation.  The 
differential  vector  connecting  the  position  r  (£a ,0  of  a  particle  with 
the  position  r  (£a  +  d£a,  c  +  ds)  of  any  neighboring  particle  is  given 
by 


dr  =  g  d?“  +  g  ds 
-a  -3 


(3.1) 


The  distance  ds  between  these  particles  is  the  magnitude  of  dr  ,  or 
ds2  =  gag  dS°  d?S  ♦  (ga3  +  g3a)  d?a  ds  ♦  g  ds2  .  (3.2) 
Initially  at  time  t  =  tQ  this  distance  is  given  as 

dS2  =  Gq6  d?“  d5S  ♦  (Ga3  +  G3a)  d?  ds  +  G33  ds2  .*  (3.3) 

A  measure  of  the  stretching  undergone  from  the  time  t  =  tQ  is 


Throughout  this  report  the  initial  values  of  variables  when  t  =  t 
will  be  denoted  by  replacing  minuscule  by  corresponding  majuscule0 
letters,  whenever  this  can  be  done  without  ambiguity.  Otherwise  a 
zero  subscript  will  be  appended  to  indicate  initial  values. 
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\  (d s2  -  dS2)  =  taB  d f  dKe  *  (Ea3  ♦  E3a)  dsa  dc  +  E33  d?2  ,  (3.4) 


where 


Eij  - 1  <hj  -  <V 


(3.5) 


This  measure  vanishes  when  no  stretching  occurs  between  two  particles. 
Moreover  when  no  stretching  occurs  in  the  neighborhood  of  a  particle, 
so  that  the  measure  of  stretching  vanishes  for  all  d?a ,  d?  ,  then  E„ 
will  vanish.  This  makes  it  reasonable  to  adopt  E. .  as  a  strain  measure. 

Due  to  the  Kirchhoff  hypothesis,  ga3  =  g3a  =  0  and  g33  =  1  for 
all  time.  Consequently,  the  strain  E_  will  assume  the  simple  form 


/  E11  E12  0 

Eij  =  (  E21  E22  0 

\  0  0  0 


(3.6) 


This  shows  that  no  stretching  or  shearing  takes  place  in  the  transverse 
direction.  Hence  we  need  only  consider  the  nonzero  components  of  the 
strain,  namely 


EoB  2  ^gaB  ~  GoB)  ' 


(3.7) 


Clearly,  the  synmetry  of  g^  implies  the  symmetiy  of  E^  . 

By  virtue  of  (2.38)  and  (2.46)^  the  nonzero  components  of  strain 
can  be  related  to  the  surface  metric  and  b  tensor: 

E„B  ■  I  t(aa6  -  AoB)  -  2C  Cba6  -  BoS)  .  (b -  B*,)]  .  (3.8) 

where  b'g  denotes  the  components  of  the  square  of  the  second  funda¬ 


mental  tensor: 


b2  =  a  .  bY  b. 


(3.9) 


From  the  definition  of  b  _  (2.15)  and  the  fact  that  n  is  a  unit  vector, 

ap 

it  follows  that  b2  can  be  written  as 
ap 

baB  =  ?’a  *  "»{S  =  '  ”  *  ?*aB  (3*10) 

Clearly  b2g  is  symmetric. 

3.2  Strain  -  Displacement  Relations 

Let  Au  be  the  displacement  undergone  by  a  middle  surface  particle 
from  some  past  time  t  to  the  current  time  t .  Denote  the  values  of 
variables  at  the  time  t-  by  an  attached  minus  superscript  and  leave  un¬ 
changed  the  values  corresponding  to  t-  Using  this  notation,  the  dis¬ 
placement  Au  becomes 


Au  =  r  -  r 


(3.11) 


The  strains  at  these  times  are 


e;6  - 1  -  v  - 2s  -  w  * <2  -  B«e>i  > 


Ea6  2  t(aaB  ■  Aoe'  '  2C  tbafi  "  BaB^  *  ^  fbafi  '  BaS^  ' 


(3.12) 


Hence,  it  follows  that  the  finite  strain  increment  is  additive: 


EaB  EaB  +  AEaB 


(3.13) 


where  the  incremental  strain  is  given  by 


4EaB  =  1  [4aaB  '  25  “aS1 


and 


AaaS  "  anB  ‘  3aS’ 


Aba&  ”  baB  baB’ 


Ab2  =  b2  -  b2: 
a6  aB  a6 


(3.14) 


(3.15) 
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In  order  to  express  the  strain  increment  in  terms  of  gradients  of 
Au  and  An  ,  we  first  observe  that 


Au  =  a  -  a 
_,a  .a  .a 


An  =  n  -  n 


(3.16) 


Hence,  from  the  definitions  of  aQg  ,  and  b2g  ,  equations  (2.6), 
(2.15)  and  (3.10),  we  have 


Aa  ,  =  a  •  Au  „  +  aQ  •  Au  +  Au  •  Au  „  , 

aC  -a  .6  _,a  ^,a  ’ 


Abo6  =  n  •  4u><j6  .  An  •  .  An  •  AUj(je  ,  (3.17) 

ale  ‘  V  ■  4n,S  *  “je  •  *;,a  *  ‘  “,B  ’ 


or  combining  differently,  we  also  have 


He  =  ?n  •  4“,s  *  ?B  ■  s?.«  '  Ha  •  He  ■ 


Ab  0  -  n  *  Au  a  +  An  ♦  a  „  -  An  *  Au  .  ,  (3.18) 

ag  _  _,ag  _  _a,6  .  ,a6  *  ^  1 

Ab2.  =  n  •  An  .  +  n  •  An  -An  •  An  „  . 

a8  _,o  _,g  _,a  >,a  _,6 


With  these  two  sets  of  expressions-  we  can  write  the  strain  increment 
in  the  following  equivalent  ways: 


He  “  2  •  He  *  ?e  ‘  H«  *  Ha  '  He1 


-2C  (n~  •  Au>ae+  An  •  aa>g  *  An  *  Au>ag) 


+C2  (n~  •  An  +  n~  •  An  +  n  •  An  )] 

-i®  -»D  ->P  -»®  -)P 


(3.19) 


and 


AE  [(a  *  Au  „  +  a0  •  Au  -  Au  •  Au  Q) 

ag  2  _a  _,B  ~8  -,a  _,a  _,8' 


-2?  (n  •  Au  aB  +  An  •  a^g  -  An  •  Au^) 


+?2  (n  *  An  +  n  •  An  -  An  „  •  An  .)] 


(3.20) 


>  _,8  “,6 


,a  ,  ,ct 


Both  these  expressions  are  exact  within  the  limits  of  the  Kirchhoff 
hypothesis.  They  are  quadratic  in  increments,  with  no  linearization 
based  on  the  smallness  of  these  increments  being  used.  An  exact  ex¬ 
pression  linear  in  the  increments  can  be  obtained  by  averaging  these 
expressions: 

1 


4E„B  =  2  l(?a  '  4“,B  *  ?8  ‘  4%> 


-2C  (n  •  Su>aB  ♦  to  •  ao>8) 


(3.21) 


with 


"?2  (",a  •  ^,8  *  ?.B  ’  4V!  • 


!a  “  7  '  9  •  "  '  1  f  >  ' 


(3.22) 


It  should  be  noted  that  the  last  quantities  have  no  significance  other 
than  as  averages,  there  not  necessarily  being  any  intermediate  position 
of  the  middle  surface  for  which  a^  forms  a  basis  and  n  a  surface  normal. 
Indeed,  ii  will  not  generally  be  a  unit  vector  nor  be  normal  to  a  . 

-  -O- 


3.3  Surface  Normal  -  Displacement  Relations 


When  the  gradient  of  the  displacement  increment  Au  is  small,  there 
is  very  little  difference  between  n~  and  n.  This  makes  the  use  of  (3.16) 
insatisfactory  for  determining  An  numerically,  since  computational  ac¬ 
curacy  could  be  impaired,  and  suggests  the  use  of  a  linear  expression 
based  on  expanding  An  as  a  series  in  Au  .  On  the  other  hand,  for 
moderately  large  Au  a  the  linear  expression  would  lose  accuracy.  Since 
our  aim  is  to  ultimately  use  this  theory  in  a  computer  code,  these  con¬ 
siderations  lead  us  to  derive  an  exact  expression  for  An  containing 

* 

Au  as  a  linear  factor. 

We  begin  by  noting  that  due  to  n  being  a  unit  vector 


An*(n  +  n)  =  (n-n)*(n  +  n)  =  0  ,  (3.23) 


and  due  to  n  being  normal  to  a^ 

An  ♦  a  =  -  n  •  a  =  -  n  •  Au  .  13.24) 

_a  _a  _  „,a 

Resolving  An  relative  to  the  basis  (a^  n),  we  obtain 

An  =  An  aa  +  Ar.  n  .  (3.25) 

a  . 

The  tangential  component. 

An  =  An  •  a  ,  (3.26) 

a  _a 

can  be  written  as  a  linear  expression  in  Au  by  using  (3.24) 

-  >ct 

An  =  -  n  ♦  Au  .  (3.27) 

a  _  „,a 

As  far  as  the  author  is  aware,  MORI. VO,  LEECH  and  WITHER  [2]  are  the 
first  to  have  addressed  themselves  to  this  problem.  The  present 
analysis  is  partially  based  on  theirs. 
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In  order  to  obtain  an  appropriate  expression  for  the  normal  componet  of 
An,  we  multiply  (3. 25)  by  (n  +  n") .  By  virtue  of  (3.23),  the  product 
vanishes  so  that 


An  n  •  a  +  (1  +  n  *  n  )  An  =  0 
a  _  ..  _ 


(3.28) 


Solving  this  equation  for  An  and  using  (3.24)  and  (3.26),  we  obtain  an 
expression  quadratic  in  An  and  hence  quadratic  in  Au  : 

Ct  •«  jCl 


aB  .  . 

a  An  An„ 

_ a  g 

1  -*  n  •  n 


(3.29) 


To  insure  that  the  denominator  does  not  get  too  small  and  thus  reduce 
coucutational  accuracy,  we  require  that  the  above  expression  for  An  be  used 
rally  when  the  mild  restriction  n  •  n  ^  0  is  satisfied,  otherwise  (3.16)2 
should  be  used. 

Rearranging  slightly,  we  can  summarize  as  follows:  when  n  •  n  >_  0  , 
determine  the  increment  An  using 


An  =  [a  + 
- 


a  1.° 

- —  n]  An 


1  +  n  •  n 


(3.30) 


,  .  .  a  aB  , 

AnQ  =  -  n  •  Au  a  ,  An  =  a  An^  ; 


(3.31) 


otherwise  when  n  •  n  <0,  use  (3.16)2-  Lastly,  had  the  basis  (aa  ,  n  ) 
been  used  to  resolve  An  ,  the  expression  for  An  would  have  assumed  the 
similar  form: 


An  =  [a~  + 
„  1  _a 


1  +  n  •  n 


where  now 


.  a 

-  n  ]  An  , 

(3.32) 

> 

0 

n 

§ 

~ 

,  0  -aB 

An  =  a  An,, 

O 

(3.33) 

I 

3 
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4.  RESULTANT  RELATIONS  FOR  MIDDLE  SURFACE  VARIABLES 


If 


Bs- 


Under  the  assumption  that  the  shell  deforms  according  to  the 
Kirchhoff  hypothesis,  we  derive  in  this  chapter  relations  connecting  the 
variables  defined  on  the  middle  surface  to  associated  variables  defined 
on  the  shell.  More  specifically,  we  obtain  the  relations  connecting 
1°  the  linear  momentum  and  spin  momentum  densities  acting  on  the  middle 
surface  to  the  particle  velocity  field,  2°  the  stress  and  couple  result¬ 
ants  to  the  stress  field,  and  S'  the  body  force  and  body  couple  resultants 
to  the  bodty  force  field  and  tractions  on  the  bounding  surfaces.  For  the 

derivation  we  also  assume  that  the  material  conposing  the  shell  is  non- 
* 

polar  and  that  the  middle  surface  variables  are  resultants  of  the  corre¬ 
sponding  quantities  defined  on  the  shell. 

4.1  Integration  Through  the  Shell  Thickness 

In  order  to  obtain  the  resultant  on  the  middle  surface  of  a  variable 
defined  on  the  shell,  we  need  to  express  the  volume  integral  as  an  inte¬ 
gral  over  the  middle  surface  of  an  integral  through  the  thickness.  We 
begin  by  letting  S  be  any  region  in  the  middle  surface  bounded  by  the 

simple  closed  curve  C.  Let  V  be  the  portion  of  the  shell  generated  by 

Yl  Q 

varying  s  between  the  limits  ±  y  as  T  ranges  over  S,  see  Figure  4.1. 

The  volume  V  is  bounded  by  two  surfaces  S  (where  £  =•  — ) ,  and  S 
h  ^  2  - 

(where  and  by  the  surface  generated  by  the  bounding  curve  C  as  5 

ranges  between  its  limits.  Using  the  basis  g^  ,  we  can  write  the  volume 

integral  of  any  variable  ^  (5°,  ?)  defined  on  V  as 


n  ** 

IIJ  i  tta>S)  dv  =  jj  [  j  2  i  g3*  dS1  d£2  ,  (4.1) 


See  NOLL  §  TRUESDELL  [10;  Sec  98]. 
*See  SGKOLNIKOFF  [7;  Eq  44.11]. 
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Figure  4.1  The  volume  V  generated  by  varying  c  between  the  limits 
±  over  the  region  S  of  the  middle  surface.  Sis  bounded  by  the 
curve  C  and  V  by  the  surfaces  Sc,  S+  and  S_. 

wherein  £2  is  the  pre-image  of  S  (see  Appendix  A).  By  virtue  of  (2.46)3 
the  last  expression  becomes 

h 

IJj  i  U“»?)  dv  =  j}  [{  2  <5  (Sa,?)  V  del  **  de1  d g2  (4.21 

1/  fi  -I 

or,  since  S  is  the  image  of  Q  , 

h 

j]J  {  dv  =  j|  [  j  2  i  (5a,e)  udg]  da  ,  (4.3) 

V  S  — 

2 

which  is  the  desired  integral  representation.  A  similar  argument  yields 
the  representation 


Figure  4.2  The  surface  S*  generated  by  varying  c  between  the  limits 

k  A  A 

±  y  over  the  arc  C.  t  is  the  unit  tangent  to  C  and  v  the  exterior  unit 

^  A  *•  ^ 

normal  to  C  tangent  to  the  middle  surface,  v*  is  normal  to  S*  and  t* 
is  tangent  to  S£  and  the  middle  surface,  but  v*  and  t*  need  not  be 
unit  vectors. 


||  £  (?“)  da  *  ||  <  aa)  y±  da* 


(4.4) 


for  any  function  £  (?a)  defined  on  either  S  or  S_  ,  with  y  and  y  the 

h  II  *r  *"  +  — 

values  of  y  at  c  =  y  and  5  =  -  y  respectively. 

We  now  seek  an  expression  relating  an  integral  over  the  bounding 
surface  S !  to  an  integral  along  the  curve  C  of  an  integral  through  the 
thickness.  We  begin  by  taking  some  connected  subset  of  C,  say  the  arc 
C,  and  the  surface  5,  generated  from  C  by  varying  ~  <  ;  <  |  ,  as  shown 
in  Figure  4.2.  Let  {  (s,  ?)  be  a  function  defined  on  S.  ,  with  s  the 

A  C 

arc  length  along  C.  The  integral  £  (s,  C)  over  S*  can  be  expressed  as 


This  expression  is  exact  only  when  h 
varying  h  it  is  a  good  approximation. 


=  const.  However,  for  slowly 


. 


(s,0  da  = 


(t  (s,c)  i  v*  I  d? 


ds  , 


(4.5) 


where  v*  is  the  (not  necessarily  unit)  outward  normal  to  given  by  the 
expression 


v*  =  t*  *  n 


(4.6) 


with 


M- .  «! 

3s  °a  ds 


5  ?o 


(4.7) 


by  virtue  of  (2.36).  It  should  be  carefully  noted  that  xa  denotes  the 
conponents  of  the  unit  tangent  t  to  C  with  respect  to  the  basis  a  .  It 
can  be  shown,  using  (2.14)^  ,  (2.46)^  and  (2.47)2  that 


?«  *"  '  -  “  ea6  ?6  ’  (4'8) 

which  when  used  in  conjimction  with  (4.7)  enables  us  to  write  v*  as 

y*  =  v  Ega  t°  gS  ,  (4.9) 


or,  in  view  of  (A. 9)  ,  as 

v*  =  y  g8  .  (4.10) 

The  reader  should  notice  that  v  represents  the  conponents  w?th  respect 

M  A 

to  a  of  the  unit  outward  normal  to  S-  along  C.  The  last  expres.ion 
«a  C 

could  be  substituted  in  (4.5),  but  we  refrain  from  this  since  no  sinpli- 
fi cation  results  at  the  moment. 

*See  S0K0LNIK0FF  [7;  p.  151]. 
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4.2  Linear  Momentum  and  Spin  Momentum  Relations 

For  the  following  derivation  we  add  to  the  Kirchhoff  hypothesis  the 
two  basic  assumptions:  1°  the  momentum  density  of  particles  composing 
the  shell  is  the  product  of  the  mass  density  and  particle  velocity  and 
2°  the  linear  momentum  and  moment  of  momentum  acting  on  any  portion  of 
the  shell  are  the  resultants  of  the  linear  momentum  and  moment  of  mo¬ 
mentum  densities  of  the  particles  composing  that  part  of  the  shell. 

We  begin  by  endowing  the  middle  surface  with  inertia  by  ascribing 
to  each  point  composing  the  surface  a  linear  momentum  P  per  unit  area 
and  a  spin  momentum  Z  per  unit  area.  Using  the  notation  of  the  last 
section,  we  have  as  a  consequence  of  applying  the  above  assumptions  to 
the  volume  V  the  two  integral  equations 

||  P  da  =  |||  p  x  dv  ,  ||  (r  *  P  *  £)  da  =  |||  x  x  (p  x)  dv  ,  (4.11) 
S  V  S  V 

with  p  the  mass  per  unit  volume  and  x  the  particle  velocity.  Upon 
applying  the  integral  representation  (4.3),  these  become 


n_ 

||  -  da  =  |j  [  |  2  P  -  V  ds]  da  » 


Jh 
S  "2 


II  6*-ll  [I, 


P  X  X  x  u 


h 

5  "  2 


(4.12) 


da  . 


Since  S  is  arbitrary,  it  follows  from  the  continuity  of  p  and  x  that  the 
integrands  are  equal: 


-lJ  pi 


x  y  d£ 


r  x  P  +  Z 


P  x  X  X  D  dc 


(4.13) 
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By  virtue  of  the  Kirchhoff  hypothesis,  (2.34)  and  (2.35),  the  last  equa¬ 
tions  become  the  momentum-velocity  relations 


p  =  Yo  V  +  Yi  1=3  ,  l  -  n  x  (Yi  V  +  Y2  u>)  ,  (4.14) 


where  Yo>  Yi»  Y2  are  the  zeroth, first  and  second  resultant  moments  of 
mass  density,  defined  by  the  general  expression 


h 

Y  =  j  2  p  v  (?)a  dc 

Jh 
"  2 


(a  =  0,  1,  2) 


(4.15) 


These  resultants  have  physical  interpretations:  yo  is  the  mass  per  unit 
middle  surface  area,  Yi  is  the  product  of  the  mass  per  unit  middle 
surface  area  and  the  distance  from  the  middle  surface  to  the  center  of 
mass,  and  Y2  is  the  moment  of  inertia  of  the  section  relative  to  the 
middle  surface.  If  the  center  of  mass  were  to  coincide  with  the  middle 
surface,  then  Yi  would  vanish  and  Y2  would  be  minimum. 


4.3  Stress  Resultant  and  Couple  Resultant  Relations 

Since  we  assume  the  material  to  be  nonpolar,  no  couple  stress  field 
is  defined  over  the  shell  and  only  a  (force)  stress  field  need  be  con¬ 
sidered.  The  force  and  torque  acting  on  any  finite  surface  element  are 
assumed  to  be  the  resultants  of  the  stress  and  the  moment  of  stress 
acting  on  the  surface  element. 


Let  C  be  any  closed  curve  in  the  middle  surface,  then  we  assume 
that  the  action  of  the  surface  exterior  to  C  on  the  part  interior  to  C 
is  equivalent  to  a  stress  resultant  vector  per  unit  arc  length 

per  unit  arc  length.  It  is  to  be 


and  a  coiale  resultant  vector 


This  is  a  statement  of  the  stress  principle  for  shells  postulated  by 
ERICKSEN  5  TRUESDELL  [11]. 
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expected  that  these  resultant  vectors  will  be  functions  not  only  of  the 
point  through  which  C  passes  but  also  of  the  curve  C  itself.  IVe  anticipate 
that  the  dependence  on  C  will  be  through  the  exterior  unit  normal  v  vo 
the  curve  tangential  to  the  surface  by  writing  the  stress  resultant  and 
couple  resultant  vectors  with  subscript  (v) . 

Using  the  notation  introduced  in  4.1  and  applying  the  above  assump¬ 
tions  to  the  surface  S~  generated  by  the  arc  C  subset  of  the  closed  curve 
C  (see  Figures  4.1  and  4.2),  we  obtain 

|  ?(j)  ds  =  ||  a  •  v  da 
C  S~ 

c  (4.16) 

|  tr  "  Q(u)  ♦  »M)  ds  =  ||  x  «  o  ■  v  da  , 

8  '  '  % 

with  o  the  stress  tensor  and  v  the  unit  exterior  normal  to  S-  ,  so  that 

-  c 

o  •  v  is  the  force  per  unit  area  over  S~  .  Since  v  is  a  unit  vector  co- 
linear  with  v*,  we  can  write 

v*  =  J  V* 1  V  ,  (4.17) 


so  that  on  applying  the  representation  (4.5)  to  the  above  integrals  they 
become 


n 

1  ?cv)  *  *  j  [  L2 °  •  r dc] 


ds 


(r x  S(v) +  45  = 


i  L 


(4.18) 


h 

2 

h 

7 


x  x  o  •  v*  d? 


ds 
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Since  these  integrals  hold  for  all  arcs  C  subsets  of  C,  it  follows  from 
the  continuity  of  the  stress  tensor  that  the  integrands  are  equal,  so 
that  in  view  of  (2.34)  we  obtain 

h  h 

Q(vj  =  |  2  a  *  v*  d£  ,  =  n  *  |  2  a  •  v*  x,  d?  .  (4.19) 

~  2  '  2 

Next  we  use  (4.10)  enabling  us  to  write  and  m^  as  linear  fimctions 
of  the  conponents  of  the  exterior  unit  normal  v  to  the  arc  C: 

h 

8(v)  =  [  Ij,2 "  ’  ?° u  dc]  u« 

2  (4.20) 

h 

?(v)  =  fc  *  jh2  -  *  g“  P  C  dC]  Va  * 

~  2 

Ct  * 

Or  defining  the  stress  resultant  tensor  Q  and  couple  resultant  tensor 
m°  as  follows 

h  h 

Q°  =  f  2  a  ♦  g“  u  d?  ,  m“  =  n  x  [  2  o  *  g“  u  £  d?  ,  (4.21) 

Jh  "  Jh  ' 

~  2  “2 

(notice  that  both  are  independent  of  the  curve  normal  v) ,  we  can  finally 

*  tt  & 

Q  and  m  are  called  tensors  because  each  is  pair  of  vectors  in  3-space 
that,  in  addition,  behave  like  components  of  vectors  in  the  (2-dimen¬ 
sional)  tangent  space  of  the  middle  surface  under  transformation  of  the 
material  coordinates  Ca.  More  properly  they  should  be  called  double 
tensor  fields,  see  [11;  Eq.  3.1]  for  a  definition. 
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t 

1 

i 

1  write  the  stress  and  couple  resultant  vectors  as 

I 

l  ?M  -  f  \  ■  ?M  ■;“'’«*  •  <4-22> 

! 

j  Because  ma  lacks  a  normal  conponent,  as  is  clear  from  (4.21)2,  we 

find  it  convenient  to  introduce  an  equivalent  tensor  Ma,  the  bending 
resultant  tensory  through  the  definition 

Ma  =  ma  x  n  ,  (4.23) 

with  the  equivalence  of  these  two  tensors  following  from  the  conjugate 
relation 


m“  =  n  x  m“  ,  (4.24) 

which  is  a  consequence  of  (4.21)2-  Like  m“,  m“  only  has  tangential 
components . 

The  stress  resultant  Q°  and  bending  resultant  M°  are  conveniently 
resolved  relative  to  the  basis  (aQ,  n)  : 

Q°  =  QaB  *  Qa  n  ,  M°  =  m“6  ag  ,  (4.25) 

and,  by  virtue  of  (4.21)  and  (4.23),  we  obtain  for  their  components 


These  relations  could  also  have  been  Stained  from  integral  equations 
of  equilibrium  by  arguments  based  on  tne  continuity  of  the  integrands 
and  the  arbitrariness  of  S,  c  f .  ERICKSEN  §  TRUESDELL  [11;  Sect.  24], 
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Writing  the  conqponents  of  the  stress  a  relative  to  the  basis  as 
follows 


lj  i  j 

o  =  g  •  a  *  gJ 


(4.27) 


we  can,  using  (2.45)^  and  (2.46)^,  replace  a  by  its  components  to  obtain 


„aB  f  2  6  ya  .  _a  [  2  3a 

Q  -J  UY«  »  *  .  Q  -]  °3  » 


(4.28) 


11 

ua8  f  2 

M  =  p 
.  h 


6  Y»  . 

!  o  11  ?  d? 


By  convention  we  call  the  membrane  components  and  Qa  the  shear 
components  of  the  stress  resultant.  The  last  relations,  as  well  as  the 
one  before,  will  hold  irrespective  of  how  the  stress  is  obtained  from 
the  strain  and/or  strain  rates,  etc.  Although  not  used  in  this  deriva¬ 
tion,  the  stress  is  assumed  to  be  symmetric 


Jt.i 1  . . . .  JWrtmillrtflUlMal 


,M  . . . . . .  t,., - - - 


4.4  External  Force  and  External  Couple  Relations 


Again,  using  the  fact  that  the  material  is  nonpolar  we  consider 
only  a  body  force  field  defined  over  the  shell.  We  assume  that  the 
force  on  a  finite  volume  is  the  resultant  of  the  stress  acting  over  the 
volume's  boundary  and  the  body  force  acting  over  the  volume,  and  that 
the  torque  on  the  volume  is  the  resultant  of  the  moment  of  stress  over 
the  boundary  and  the  moment  of  body  force  over  the  volume. 

We  begin  by  accounting  for  the  influence  of  the  external  environ¬ 
ment  on  the  middle  surface  by  assigning  to  each  point  on  the  surface  an 
external  force  F  per  unit  area  and  an  external  couple  c  per  unit  area. 
We  then  apply  the  above  assumptions  to  the  volume  V  bounded  by  the 
surfaces  S  ,  S+  and  S  (see  Figure  4.1)  to  obtain,  after  taking  care  to 
eliminate  the  contributions  from  by  applying  (4.16)  to  C,  the  two 
integral  equations 


||  F  da  =  HI  b  dv  +  ||  t+  da  +  ||  t  da  , 


S  S 

+ 


(4.30) 


||  (r  x  F  +  c)  da  =  |||  x  x  b  dv  +  ||  x+  x  t+  da  +  ||  x  *  t  da  , 


with  b  the  body  force,  and  t  and  t  the  tractions  on  the  surfaces  where 

-*+  -  ** 

C  =  y  and  £  =  -  2  ,  respectively.  Since  the  exterior  unit  normals  to 
S+  and  S_  are  n  and  -n,  we  obtain  for  the  tractions 


t+  =  a  •  n  ,  t '  =  -  a  ♦  n 


(4.31) 


When  we  apply  (4.5)  and  (4.4)  to  the  last  integrals,  they  become 


49 


50 


*  L2  “8  bB  U  dC  *  “*8  '♦  "*  *  v-l  l- 

•  n 


F  =  [  2  b3  v  d;  +  t|  u+  +  t3  \i_ 
Jh 
'  2 


(4.36) 


:a  =  j  2  b3  V  C  dc  +  |  (y+“  t3  y+  -  u_“  t3  yj  , 


wherein  the  conponents  of  b  and  t  are  with  respect  to  the  basis  gi  : 


b1  -  g1  •  b  ,  t1  =  g1  •  t 


(4.37) 


'  I  Ml 


5.  EQUATIONS  OF  MOTION  OF  A  SHELL 


In  this  chapter  we  first  derive  general  equations  of  notion  for  a 
shell,  with  no  assumptions  (such  as  the  Kirchhoff  hypothesis)  being  made 
about  the  deformation  through  the  thickness.  Next  the  consequences  of 
mass  being  conserved  on  the  resultant  moments  of  mass  density  are  found. 
These  are  then  conbined  with  the  resultant  relations  of  the  last  chapter 
to  simplify  the  equations  of  motion.  Lastly,  the  SANDERS-LEONARD  proce¬ 
dure  is  used  to  eliminate  the  shear  components  of  the  stress  resultant, 
resulting  in  a  single  vector  equation  of  motion. 

5.1  Derivation  of  Equations  of  Motion 

The  derivation  is  based  on  formulating  the  principles  of  balance  of 

linear  momentum  and  moment  of  momentum  for  the  reference  surface  in  inte- 

★ 

gral  form.  From  the  last  chapter  we  have  that  the  reference  surface  is 

endowed  with  linear  momentum  P  and  spin  momentum  l  per  unit  surface  area. 

Also  the  mechanical  influence  of  the  external  environment  on  the  reference 

surface  is  specified  through  the  external  force  F  and  external  couple  c 
**  - 

per  unit  area. 

Lastly,  we  have  from  the  previous  chapter  that  for  any  simple  closed 
curve  in  the  reference  surface,  the  mechanical  influence  of  the  part  of 
the  surface  exterior  to  the  curve  on  the  part  interior  is  transmitted 
through  the  stress  resultant  vector  Q^-j  and  the  couple  resultant  vector 
m^  per  unit  arc  length.  Using  these  variables,  we  can  state  the  princi¬ 
ples  of  balance  of  linear  momentum  and  moment  of  momentum  for  any  simply 
connected  region  S  of  the  reference  surface  with  simple  closed  boundary 
C  (see  Figure  5.1)  in  the  following  integral  forms: 


We  use  reference  surface,  rather  than  middle  surface,  to  emphasize  that 
the  results  of  this  section  are  independent  of  the  Kirchhoff  hypothesis. 

For  completeness,  we  should  introduce  here  the  contact  force  Q  and  con¬ 
tact  couple  m  per  unit  arc  length  to  account  for  the  influence  of  the 
environment  entering  through  the  physical  boundaries,  but  we  postpone 
this  until  the  next  chapter,  since  these  quantities  are  not  used  here. 


PRECEDING  PAGE  BUNK 
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Figure  5.1  The  action  on  the  region  S  of  the  external  force  F, 
external  couple  c,  stress  resultant  vector  and  couple  result 
vector  m,  x.  F  and  c  are  shown  acting  on  c  typical  point  p  in  S 


vector  m(vj-  F  and  c  are  shown  acting  on  c  typical  point  p  in  S 
located  by  position  vector  r,  while  and  are  shown  acting  on 
a  typical  point  on  the  boundary  C  of  tfie  surface  S.  Also  indicated  is 
the  relative  orientation  along  C  of  the  surface  unit  normal  n,  the 
boundary  unit  tangent  t  and  external  unit  normal  v. 


. 


at  If !  *■  *  f  Sm  *  *  JJ  l da  • 


(5.1) 


||  (r  x  p  +  l)  da  =  (I  (r  *  Q,v-  +  m^)  ds  +  ||  (r  *  F  +  c)  da  . 


From  these  integral  relations  we  derive  the  differential  equations 

of  motion.  Assuming  that  the  inertia  terms  in  the  left  hand  integrals 

* 

are  continuous,  we  can  apply  Reynold's  transport  theorem  to  bring  the 
material  derivative  past  the  integral  signs: 


11 


+  V*v  P)  da  =  0  ds  +  ||  F  da 


7  ?(v)  JJ 

C  ~  S 


J j  l  +  V-v  l  +  v  x  P  +  r  x  CP  +  V-v  P)  da  =  (5.2) 

'  '  i  J  . 


r  *  8(v)  +  ®(v)5  ds  +  ||  (r  x  F  +  c)  da  . 


As  noted  earlier,  equations  (4.22)  relating  and  m^  to  tensors 
9°  and  ma  can  be  obtained  directly  from  the  balance  equations  (5.1) 
without  recourse  to  the  special  material  and  geometric  assumptions  of 
the  last  chapter.  Assuming  only  that  the  surface  integrands  are  bounded 
and  the  line  integrands  continuous,  the  usual  limit  argument  using  a 


A  derivation  of  Reynold's  transport  theorem  for  a  surface  is  included 
in  Appendix  A. 

k 

See  footnote  on  page  47. 
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sequence  of  regions  in  S  gives  (4.22),  but  without  the  restriction  that 

ma  be  tangential.  Substitution  of  (4.21)  in  the  above  integrals  allows 
—  * 

us  to  use  the  divergence  theorem  to  replace  the  line  integrals  by 
surface  integrals: 


|j  (P  +  V-y  P)  da  =  JJ  (Qa, q  +  F)  da 
S  S 

j|[£.V.vf.»^.rv(P.v.v  P) 

s 


da 


(5.3) 


11  [*. 


a  x  Qa  +  c  +  r  *  (Qa.  +  F) I  da  . 

.a  ,  _  _  _  I  a  1 


:>J 


Finally,  assuming  the  above  integrands  continuous,  then  by  virtue  of  5 
being  arbitrary  it  follows  that  (5.3)  can  be  satisfied  if  and  only  if 


P  +  V*v  P  =  Q  ,  +  F 

_  _  '  la 


£  +  V*v  £  +  v  *  P  =  ma.  +  a  *  Qa  +  c 
__  _  _  _  I  a  „a  _ 


(5.4) 


These  are  the  differential  equations  of  motion  for  a  shell,  the 
first  enfcodying  the  principle  of  balance  of  linear  momentum  and  the 
second  the  principle  of  balance  of  moment  of  momentum.  Under  the 
assumed  continuity  conditions,  they  are  fully  equivalent  to  (5.1).  They 
are  independent  of  the  properties  of  the  material  conposing  the  shell 


A  derivation  of  the  divergence  theorem  for  a  surface  is  included  in 
Appendix  A. 
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and  of  the  Kirchhoff  hypothesis,  with  no  restrictions  placed  on  m  and 
c  being  tangential.  Setting  the  inertia  terms  on  the  left  hand  sides 
equal  to  zero,  we  get  the  differential  equations  of  equilibrium  for  a 
shell: 


Q  .  +  F  =  0  , 

1  la 


cl  _a 

m.  +a  k  Q  +  c  =  0 
_  1  a  _a  j 


(5.5) 


They  are  the  vectorial  forms  of  the  equations  derived  by  ERICKSEN  § 
TRUESDELL  [11;  Eq.  25.2.,. 

Next  we  write  the  equations  of  motion  in  terms  of  components  normal 
and  tangential  to  the  reference  surface.  We  begin  by  resolving  the 
variables  in  (5.4)  into  normal  and  tangential  components: 


P»PPa6.Pn  ,  t  =  , 


na  „aB  .  na  a  aB  a 

Q  =  Q  ag  +  Qn,  m  =  m  afi  +  mn  , 


F  =  F6  ae  *  F  n  ,  c  =  c8  a6  .  c  a  . 


By  virtue  of  (2.26),  (2.29)  and  (2.33),  it  follows  that 


(5.6) 


P  =  [P£  +  PY  (ve(  -  bB  v)  +  P  u>B)  ag  +  [P  -  P6  M  )  n  , 


As  a  consequence  of  (2.20),  we  have 


8V  tQaVb!  Qal  +  [Q“t  a  +  baB  Qa6ln 


(5-8) 


a  ,  aB  ,8  a.  ,  a  ,  aB, 

?  la  =  tm  la  '  ba  m  ]  ?B  +  [m  la  +  baB  m  ]  " 


The  cross  product  terms  in  (5.4)  are  resolved  by  using  (2.14): 


v  * p  ’  <■*  p“  -  v“  p>  ?B  *  £»6  v° pB  "  ■ 


(S.9) 


a  x  Qa  =  -  e  Qa  a8  +  e  Qa8  n 
~a  aB  -  aB  . 


Substitution  of  the  last  three  results  and  (2.27)  into  the  equations  of 
motion  (5.4)  yields 


P8  +  Pa  (v8  -  b8  v)  +  P8  (v“  -  b°  v)  +  P  io8  =  Qa8,  -  bB  Qa  +  F8  , 

la  a  v  I  a  a  J  x  ~  ’ 


I  a  a 


P  +  P  (va.  -  ba  v)  -  PQ  u)  =  Qa.  +  b  „  Qa8  +  F  , 
la  a-’  a  x  I  c*  aB  H  * 


£6  +  ta  (v8  -  b8  v)  +  £8  (va.  -  ba  v)  +  l  a)8  +  e  8  (v  Pa  -  v“  P) 


(5. 


I  u 

I  a  a 


1  w 

I  a  a 


aB  ,  S  a  B  na  ^  B 
=  m  .  -bm-e  Q+c 
la  a  a  ^ 


5  „  ,  a  .  a  .  »a  a  nB  a  ,  aB  »aB 

t  +  t  (v  ,  -b  v)  -  £  oj  +  e.v  P  =  m  .  +bflm  +  e  „  Q  +  c 
v  I  a  a  '  a  aB  I  a  aB  aB 


the  tangential  and  normal  components  of  the  equations  of  motion  for  a 


s 

l 

\  shell.  The  first  and  second  equations  embody  the  principle  of  balance 

of  linear  momentum  and  the  third  and  fourth  the  principle  of  balance  of 
moment  of  momentum.  On  setting  the  inertia  terms  equal  to  zero,  the 
component  form  of  the  equations  of  equilibrium  result,  coinciding  ex¬ 
actly  with  those  derived  by  ER1CKSEN  5  TRUESDELL  [11;  Eq.  (26.6  -  26.9)]. 

An  alternative  formulation  of  the  equations  of  motion  (5.4),  in 
which  the  covariant  derivatives  are  replaced  by  partial  derivatives  and 
the  velocity  divergence  terms  vanish,  can  be  obtained  by  utilizing  the 
relations  (2.18)  and  (2.28).  Defining  the  new  variables: 

P*  =  it2  P  , 

l*  =  i  , 

it  follows  by  virtue  of  these  relations  that  the  equations  of  motion 
become: 

Q*“  +  F*  ,  l*  *  v  x  P*  =  m*“  +  a  x  Q*a  +  c*  .  (5.12) 

_  » a  __  _  ’a  _a  _ 

All  equations  derived  in  this  section, as  mentioned  before,  are 
completely  general  and  do  not  depend  on  the  form  assumed  for  the  defor¬ 
mation  through  the  shell  thickness.  Hence,  these  equations  are  valid 
for  all  shell  theories. 

5.2  Conservation  of  Mass 

We  now  derive  the  consequences  of  mass  being  conserved  on  the 
resultant  moments  of  mass  density  Yo»  Yi>  Y2'  As  a  result  of  the  mass 
being  conserved  for  any  portion  of  the  shell',  the  mass  density  p  satis¬ 
fies  the  well  known  equation  of  mass  continuity 

p  +  p  Div  x  =  0  .  (5.13) 


Q*“  -  aH  Qa  , 


.  U  ^2  U 

m*  =  a  m 


% 

F*  =  a2  F 


c*  =  a*5  c 


(5.11) 
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Writing  the  divergence  in  the  (£°,£)  coordinate  system: 

8x  * 

Div  x  =  g“  •  x,a  +  g3  ♦  ~  ,  (5.14) 

we  obtain  by  virtue  of  (2.29),  (2.36)  and  (2.45) 3 

Div  x  =  g°  *  ga  +  n  •  tu  .  (5.15) 

Since  w  is  tangential,  the  second  right  hand  term  vanishes  and  using 
(2. 45) 2  we  have 

.  •  1  ct8  •  k  ,  Jc 

Div  x  =  j  g  P  gag  =  g2  /  g2  •  (5.16) 

Applying  (2.46) 3  to  this  result,  the  equation  of  mass  continuity 
becomes 

ft  (P  u  a*'2)  =  0  .  (5.17) 

Consequently,  the  product  p  y  a2  is  a  time  constant  and  hence  equal  to 

its  initial  value  at  time  t  =  t  : 

0 

h  k  ** 

p  p  a2  =  Po  yo  A2  .  (5.18) 


See  GREEN  5  ZERNA  [12;  Eq.  1.12.36]. 
See  the  footnote  on  page  3i. 
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Applying  this  result  to  the  defining  equations  for  Yfl  (4.15),  gives  us 


u  i' 

a^a  =  A2la 


(a  =  0,1,2) 


(5.19) 


with  r  the  initial  values  of  v  .  Also,  the  new  variables 
a  a. 


v*  =  a  y 
T  a  'a 


(a  =  0,1,2) 


(5.20) 


are  time  constants  and  their  material  derivatives  must  vanish: 


y*a  =  0 


(a  =  0,1,2) 


(5.21) 


From  (2.28)  we  have  that  the  last  result  is  equivalent  to 

Ya  +  V‘!Ya*°  (a  =0,1,2)  .  (5.22) 


We  can  now  express  the  general  inertia  terms  in  the  equations  of 

•  •.  •  *• 

motion  (5.4)  in  terms  of  the  two  accelerations  v  =  x  and  w  =  n  ;  for 
conbining  the  last  result  with  the  momentum  relation  (4.14),  we  obtain 


P  +  V*v  P  =  Yn  v  +  Yi  U 


l  +  V*v  l  +  V  x  p  =  n  x  (yi  V  +  Y2  “) 


(5.23) 


With  these  equations,  the  equations  of  motion  can  be  written  as: 
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(5.24) 


Y0  v  +  Yi  u  =  Qa(o  +  F 


n  X  (Yl  v  +  Y2  “)  =  m 


I  a 


_a 

a  x  Q  +  c 
.a  _ 


The  form  assumed  by  the  inertia  terms  in  the  above  equations  is  a 
consequence  of  the  Kirchhoff  hypothesis.  Another  consequence  of  the 
Kirchhoff  hypothesis  is  that  m°  and  c  must  be  tangential,  as  indicated 
by  (4.21)2  and  (4. 33) 2-  Taking  this  result  into  account  by  replacing 
m°  and  c  by  the  equivalent  M°  and  C,  the  second  of  the  above  equations 
becomes 


n  x  (Yl  v  +  y2  «)  =  n  x  ag  Cm“B,  a  +  CB  -  QB) 

(5.25) 

♦  aa  x  ag  (QaB  -  b“  MyB) 

This  equation  is  easily  separated  into  conponents ,  giving  for  the 
normal  conponent 


0  =  e  (QaB  -  b°  MyB) 
aB  Y 


(5.26) 


and  for  the  tangential  components 


aB  •  (yi  v  +  Y2  u)  =  Ma6|0  +  CB  -  QB 


(3.27) 


These  equations  express  the  balance  of  moment  of  momentum  in  the  normal 
and  tangential  directions,  respectively.  From  (4.28)1  3  using  (2.38), 
we  have  the  identity 
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(5.28) 


h 

<f6  -  b°  =  f  2  U*  uj  a"6  u  d, 
y  r  Jh  Y  6 

"  2 

which  shows  that  if  oaS  is  symmetric  then  (S.26)  is  identically  satis¬ 
fied.  Hence,  the  tangential  components  of  the  stress  tensor  being 
symmetric  implies  the  balance  of  moment  of  momentum  in  the  normal  di¬ 
rection,  making  (5.26)  redundant.  For  this  reason,  (5.26)  is  often  not 
included  among  the  equations  of  motion  of  Kirchhoff  shells.  In  the 
next  section  we  shall  show  how  this  equation  can  be  utilized  to  affect 
a  symnetrization  of  the  membrane  and  bending  resultants. 

The  consequences  of  mass  conservation  can  also  be  applied  to  the 
equations  of  motion  for  the  starred  variables  (5.12).  Combining  (4.14), 
(5.11)  and  (5.21),  we  obtain 

P*  =  y*0  v  +  y*1  u  ,  +  v  *  P*  =  n  x  (Y*i  v  +  Y*2  m)  (5.29) 


Using  the  last  expressions  to  replace  the  inertia  terms  in  (5.12),  we 
have 


Y*0  V  +  Y*i  w 


+  F* 


(5.30) 


n  X  (y*l  v  +  Y*2  <*0  =  », 


a  x 
_a 


+  c* 


Replacing  m*a  and  c*  in  the  second  of  the  above  equations  by  the  equi¬ 
valent  tangential  vectors 


*See  for  example  KRAUS  (13;  Eq.  (2.71)]. 


this  equation  becomes 


n  x  (Y*!  v  +  y*2  «)  =  n  x  ag  (M*aS,rt  +  I^8  M*ay  *  C*8  -  Q*8) 


a  ay 

>*a6  _  m*y8. 


(5.32) 


+  a  x  ac  (Q*  -  b 

~oc  Y 


giving  for  its  normal  component 


m 

WE 


°  =  ea6  (Q*a3  '  by  M*Y3)  (5-33) 

and  for  its  tangential  components 

a8  •  (y*i  v  +  y*2  «)  =  M*a6,  +  T  8  M*“Y  +  C*8  -  Q*8  .  (5.54) 

**  ■*  —  o  u  Y 

Equations  (5.29),  (5.30),  (5.32)-(5.34)  in  the  starred  variables 
correspond  to  equations  (5.23)-(5.27)  in  the  unstarred  variables. 

Again,  (5.33)  is  implied  by  the  symmetry  of  the  tangential  components 
of  the  stress  tensor  and  in  that  sense  is  redundant. 

The  various  forms  of  the  equations  of  mass  conservation  and  of  the 
equations  of  motions  obtained  in  this  section  are  based  on  the  premise 
that  the  shell  deforms  according  to  the  Kirchhoff  hypothesis.  This 
does  not  mean  that  the  resultant  relations  obtained  in  Chapter  4  for 
Kirchhoff  shells  need  necessarily  be  used  with  these  equations,  but  only 
that  the  equations  are  compatible  with  these  relations.  Other  resultant 
relations  can,  of  course,  be  used.  However,  it  is  clear  from  the  special 
forms  of  the  inertia  terms  and  the  tanger.cy  of  the  couple  resultant  ma 
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or  m*a  and  the  external  couple  c  or  c*  to  the  middle  surface  that  the 
equations  of  motion  derived  in  this  section  cannot  be  as  general  as 
those  obtained  in  section  5.1. 

5.3  Reduced  Equation  of  Motion 

Whenever  the  Kirchhoff  hypothesis  is  used,  it  is  necessary  to  dis¬ 
regard  the  values  given  by  the  resultant  relations  (4.28).  _  for  some 
components  of  the  stress  resultant  Q  .  The  reason  for  this  is  that  with 
the  Kirchhoff  hypothesis  the  use  of  all  the  equations  of  motion  and  all 
the  resultant  relations  yields  an  overdeterminate  theory.  To  make  this 
clear,  assume  that  at  some  instant  the  position  r,  the  velocity  field 
v,  the  stress  field  a,  body  force  field  b  and  traction  fields  t+  6  t 
are  known.  Through  relations  (4.21)  §  (4.33)  the  resultants  Q®,  ma,  F 
§  c  will  be  known  and  consequently  the  right  hand  side  of  the  equations 
of  motion  (4.21)  completely  determined.  As  a  consequence  of  condition 
(2.31)  making  u>  a  function  of  V|  ,  the  equations  of  motion  become  a  set 
of  five  scalar  partial  differential  equations  for  the  three  components 
of  the  acceleration  v  .  Since  the  five  equations  are  independent,  this 
clearly  implies  that  v  is  overspecified. 

It  is  usual  to  disregard  the  values  of  the  shear  components  Qa  given 
by  resultant  relations  (4.28)2  and  take  as  significant  the  values  satis¬ 
fying  the  tangential  components  of  the  moment  of  momentum  equation  (5.27). 

This  is  accomplished  by  using  (5.27)  to  eliminate  the  shear  components 

* 

from  the  linear  momentum  equation  (5.24) Since,  as  mentioned  earlier, 
(5.26)  is  redundant  and  can  be  disregarded,  with  the  elimination  of  Q 
the  equations  of  motion  are  reduced  to  a  single  vector  equations  involv¬ 
ing  the  components  and  M®*5. 


KRAUS  [13;  P.47]  indicates  the  use  of  this  elimination  without  pointing 
out  its  necessity.  MORINO,  LEECH  and  WITHER  [2;  P.19]  do  mention  the 
need  for  this  elimination  and  go  on  to  use  it  in  developing  their 
equations. 
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Rather  than  following  this  procedure  to  reduce  the  equations  of 
motion,  we  shall  use  a  more  elegant  method  and  use  both  (5.26)  and  (5.28) 
to  eliminate  the  components  Q°  and  (/  from  the  linear  momentum 
equation  (5.24)^.  In  doing  this,  we  follow  SANDERS  [14]  and  LEONARD  [15] 
by  defining  modified  membrane  components,  but  differ  from  these  authors 
in  that  we  employ  the  equations  of  motion  rather  than  the  equations  of 
equilibrium  and  we  obtain  a  vectorial  rather  than  a  componental 
representation . 

We  begin  by  defining  the  modified  membrane  components  of  the  stress 
resultant: 


QaB  =  Q°6  +  b6  MYa  =  Q6“ 


(5.35) 


wherein  symmetry  follows  from  (5.26).  With  the  aid  of  (2.20)3>  we  use 
(5.27)  and  (5.35)  to  express  the  stress  resultant  as  follows: 


Q°  =  (MBa  n), p  +  n“e  ag  +  Ca  n  -  a“  •  (Yi  v  +  yz  w)  n.  (5.36) 


Substituting  this  expression  into  the  linear  momentum  equation  (5.24)^ 
and  using  the  identity 


(5.37) 


We  adapt  the  practice  of  enclosing  a  pair  of  indices  in  a  bracket  or 
parenthesis  to  denote  the  antisymmetric  or  symmetric  part  with  respect 
to  these  indices. 

ft 

This  identity  easily  follows  from  the  observation  that  the  anti¬ 
symmetric  part  of  any  two  index  quantity  ^aB  can  be  written  as 

t*1  -  ^  <•  »  <"?>  ■  <iaS  ■  <oS  «,oB  -  r*  '  )  .  o . 
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we  obtain  the  reduced  equation  of  motion 


Y0  V  +  Yi  to  +  [aa  •  (Yi  V  +  Y2  u)  n]la 

•  CMB“  ?)|6a  >  M3S  h  *  c“  ”),«  *  l 

with  MaB  denoting  the  symmetric  part  of  the  bending  components: 
M°6  =  M<°»  =  i  (m“6  *  MBa)  . 

The  abbreviations 

=  =  Q“6  aB  ,  C“=C*» 

permit  us  to  write  the  reduced  equation  of  motion  more  compactly: 
Y0  v  +  Yi  u  +  [a“  •  (Yi  V  +  Y2  <*>)  n](a 

-  ^Vsa  +  sX  +  ?  • 

Setting  the  inertia  terms  on  the  left  hand  side  equal  to  zero,  we 
the  reduced  equation  of  equilibrium 

which  is  the  vectorial  representation  of  the  equations  derived  by 
SANDERS  [14,  Eq.  55  $  56].  ■ 


(5.38) 


(5.39) 


(5.40) 


(5.41) 


obtain 


(5.42) 
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The  elimination  just  performed  reduces  the  equations  of  motion  from 
two  vector  equations  (5.24)  (six  scalar  components)  to  a  single  vector 
equation  (5.41) (three  scalar  components) .  Simultaneously,  the  ten  vari¬ 
ables  Qa^,Qa  §  ma^  have  been  replaced  by  six  variables  Qa^  §  M°^.  Had 
only  the  shear  components  been  eliminated  as  indicated  before ,the  thus  re¬ 
duced  equation  would  have  in  wived  the  eight  variables  Qag  and  Ma^.  We 
will  discuss  this  further  in  section  8.1. 

',ng  '"‘aft 

The  resultant  relations  for  the  values  of  Q  and  M  are  easily 
determined  by  combining  (5.35)  and  (5.39)  with  (4.28^  ^  >  giving  us 


f  2  r  a  yB  B  yo 

Jh  [V 


2 


a  B  Y^-.  j 
UY  u6  o  ]  v  dc 


h 

2 

h 

2 


[*•;  o 


YB 


B  Y  a, 
vr  crT  ] 


u  C  dc 


(5.43) 


We  note  that  the  symmetry  of  Qa^  follows  from  the  symmetry  of  aa^  , 
which  is  not  surprising  in  view  of  the  redundancy  of  (5.26)  due  to  (5.28). 

We  conclude  by  reformulating  the  reduced  equation  of  motion  using 
the  starred  variables: 


Q*a  =  a*5  Qa  ,  M*°6  =  a2  Ma6  ,  C*a  =  a2  Ca 


(5.44) 


With  the  aid  of  (2.18)  we  can  express  (5.41)  as  follows 


Y*o  v  +  y*i  u  +  [a  •  (y*i  v  +  y*2  u>)  n] . 


=  +  N*a],o  +  F* 


(5.45) 
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In  this  chapter  we  obtain  natural  boundary  conditions  for  the  re¬ 
duced  equation  of  motion.  We  do  this  through  a  variational  procedure 
based  on  the  principle  of  virtual  work.  The  variation  in  the  rotation 
of  a  shell  element  is  assumed  to  satisfy  the  Kirchhoff  hypothesis.  The 
three  classic  boundary  conditions  of  the  free,  hinged  and  clamped  edge 
are  obtained  along  general  boundary  curves  and  then  specialized  to  the 
coordinate  curves.  The  boundary  conditions  along  coordinate  curves  for 
a  symmetry  edge  are  also  obtained. 

6.1  Variation  of  Strain  Energy 

Our  goal  is  to  find  the  natural  boundary  conditions  associated  with 
the  reduced  equation  of  motion,  using  the  principle  of  virtual  work. 
However,  since  our  concern  is  with  deriving  boundary  conditions  rather 
than  initial  conditions ,  we  can  disregard  the  inertia  terms  and  work 
just  as  well  with  the  reduced  equation  of  equilibrium  (5.42).  Therefore 
we  can  state  our  objective  more  precisely  as  follows:  to  formulate  a 


tea  i  nWnItai 


contact  couple  per  unit  arc  length  applied  along  the  boundary  C  and,  as 
already  mentioned,  F  and  c  are  the  external  rorce  and  external  couple 
per  unit  area  over  the  surface  S. 

The  rotational  part  60  of  the  variation  in  the  basis  (aQ,  n)  is 
given  by  the  expression 


60  =  •-  (aa  x  6a  +  n  *  6n) 
~  2  -  _a  J 


(6.2) 


This  expression  assumes  no  dependence  of  6n  on  da^,  even  though  n  is 
normal  to  aQ.  We  set  up  such  a  dependence  by  requiring  that  the  vari¬ 
ation  in(aQ,  n)  satisfy  the  Kirchhoff  hypothesis.  Consequently  the 
variations  6n  and  6a^  must  be  such  that  the  magnitude  of  n  is  unchanged 


This  expression  is  easily  derived  by  considering  any  vector  ^  with 
fixed  components  relative  to  the  basis  (aa,  n) : 


^=^aaa  +  ^n  =  ^‘aCtaa+^’nn 


Its  variation  due  to  varying  (a^,  n)  is 


6((  =  (5°  6aa  +  |J  on  h  •  aa  da^  +  $  •  n  5n 


and  the  rotational  part  of  6 ^  is  given  by  the  antisymmetric  part: 


60  x  1$  =  i  i  •  (a“  6a  -  5a  aa  +  n  6n  -  6n  n)  , 
^  _  2  _  _  _a  .a  .  .  _  '  ’ 


from  which  the  above  expression  follows. 


is^Sslvi 
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and  n  remains  normal  to  \  during  the  variation,  which  gives  us  the 
conditions 

n  •  6n  =  0  ,  a  •  6n  +  n  •  6a  =  0  .  (6.3) 

.  _  _a 

Combining  these  conditions  with  the  expression  for  66,  we  have  as  a 
consequence  of  the  Kirchhoff  hypothesis  that  66  is  completely  determined 
by  6a  : 

66  =  4-  aa  x  (6a  +  n  n  •  6a  )  .  (6.4) 

Our  aim  is  to  find  an  appropriate  expression  for  6V;  one  which  when 
substituted  in  the  principle  of  virtual  work  (6.1)  yields  the  reduced 
equation  of  equilibrium.  We  begin  with  the  requirement  that  the  princi¬ 
ple  of  virtual  work  be  consistent  with  the  vanishing  of  the  expression 

«A  =  (Q“|a  +  F)  •  6r  ♦  (m“  a  ♦  aa  x  Qa  +  c)  .  66  .  (6.5) 

Motivation  for  this  requirement  follows  from  the  fact  that  <5 A  does 
indeed  vanish  when  the  equations  of  equilibrium  closest  to  physical 
reality  (5.5)  axe  satisfied.  Our  analysis,  however,  docs  not  depend 
upon  (5.5)  being  satisfied,  but  only  on  (6.5)  vanishing. 

A  change  in  the  order  of  differentiation  in  the  last  expression 

gives 

6A  =  (Qa  «  6r  +  ma  •  68),  +  F  •  6r  +  c.  •  66 

(6.6) 

-  (Qa  •  6aa  -  aa  x  Qa  .  66  +  ma  •  56,  )  . 
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With  the  aid  of  (2.16)  and  (2.20),  the  derivative  of  69  can  be  written 
as  follows 


69.  =  4  [n  x  ae  (bY  6a  .  -  2  6b  _)  +  a6  *  a  6r  .Y]  ,  (6.7) 

,la  2  ,  a  y6  aB  ,  _y  <*B 


which  when  combined  with  (4.24)  gives  us 


a  1  .a  wyB  r  wa8  r, 

m  *  69.  =  —  b  M  6a  .  -  M  6b  . 

a  2  y  aB  a6 


(6.8) 


We  also  have  from  (6.4)  the  identity 


Qa  •  6a  -  a  x  Qa  .  69  =  i  QaB  6a  c 
,  ,a  .a  _  2  aB 


(6.9) 


In  view  of  definitions  (5.35)  and  (5.39),  substitution  of  the  last  two 
expressions  in  (6.6)  results  in 


6A  =  (Qa  •  6r  +  ma  •  69) j q  +  F  •  6r  +  c  •  69 


-  d  QaB  6a  .  -  MaB  6b  ) 
2  ats  aB 


(6.10) 


Integrating  this  equation  over  the  surface  S  and  applying  the  divergence 

* 

theorem  for  surfaces,  we  obtain 


’  r 

6A  da  =  0  (Q  •  6r  +  ma  *  66)  vq  ds  +  J  (F  •  6r  +  c  •  69)  da 


-  ff  4 

JJ 


(6.11) 


QaB  6a  ,  -  MaB  6b  e)  da  . 
aB  aB 


See  Appendix  A  for  a  derivation  of  this  theorem. 
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Comparison  of  this  result  with  (6.1)  shows  that  with 


5V  =  y  QaS  6a  -  Ma6  6b 
2  x  afi  a  8 


(6.12) 


the  principle  of  virtual  work  and  the  vanishing  of  6A  are  consistent. 

We  have  yet  to  show  that  this  choice  of  6V  implies  the  reduced 
equation  of  equilibrium  (5.42),  although  the  appearance  in  (6.12)  of 
precisely  those  mechanical  variables  occurring  in  the  reduced  equation 
strongly  suggests  such  a  connection.  Returning  to  (6.5),  we  use  (4.24) 
and  (4.34)  to  replace  ma  and  c  by  Ma  and  C,  giving  us 


6A  =  [Q°\  +  F]  •  6r  +  [n  x  (Ma.  +  C)  +  a  x  (Qa  -  b“  MB)]  •  66  .  (6.13) 

*  I  u  •»  —  —  —  l  (Z  _  Ot  „  P-.  ~ 


Using  (6.4)  to  eliminate  60  and  rearranging,  this  becomes 


=  [Q2,  a  *  F]  •  6r  -  [(MB°  n),  p  +  (Q(aB)  ♦  b<“  MyB))  ag 


+  C°  n  -  Qa]  •  6a 


(6.14) 


Using  the  fact  that  6aQ  =  6r|  q  to  change  the  order  of  differentiation, 
we  obtain,  in  view  of  definitions  (5.35),  (5.39)  and  (5.40),  the  result 


6A  =  [(M8cu  +  Q-1  +  C°)|o  ♦  F]  •  6r 


-  n)|B  v  Q°  +  c“  -  Qa]  .  6r}j 


Confining  (4.24)  and  (6.4)  tc  obtain 


ma  •  60  =  -  M°B  n  •  6a. 
'  -  -p 


(6.15) 


(6.16) 


and  substituting  this  result  and  the  previous  expression  for  6A  in 
(6.11),  we  obtain  after  rearranging  terms  the  integral  identity 


which  on  identifying  6V  as  per  (6.12)  becomes 

||  SV  da  +  ||  [(MSa,B  +  Q°  -  C“)|a  ♦  F]  •  6r  da 
5  S 

=  |  (F  •  6r  +  c  *  60)  da  +  <j>  {[(MB“  n))g  +  Q°  +  C°]  *  6r  (6.18) 

5  C 

-  M°8  n  *  6a  }  v  ds 
_  _p  u 

From  comparing  this  result  with  (6.1),  we  conclude  that  the  principle  of 
virtual  work  implies  the  integral  equation 

[  [  (MSol(  b  +  Qa  +  Ca)la  +  F]  •  6r  da  +  0  (Q  •  6r  +  m  •  69)  ds 
s  C  (6-19) 

=  |  {[(M6a  n),  g  -  Qa  -  C°J  *  -  M°8  n  *  oa,}  va  ds. 

r> 
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Since  6r  is  arbitrary,  we  first  set  6r  =  0  over  S,  giving  us 


(Q  •  6r  +  m  •  66)  ds  =  <j>  {[(MP  n),  g  +  Q  +  C  ]  •  6r 


(6.20) 


-  M  D  n  •  6a,}  v  ds 
-C>  a 


It  then  follows  from  (6.19)  that  for  arbitrary  or 


||  [(MB“1$  +  Q“  +  C“)la  +F]  •  6r  da  =  0 


(6.21) 


Assuming  that  the  integrand  of  this  expression  is  continuous,  the  arbi¬ 
trariness  of  6r  implies  that  the  integrand  itself  must  vanish,  which  is 
clearly  equivalent  to  the  reduced  equation  of  equilibrium  (5.42)  being 
satisfied.  Hence  we  have  shown  that  with  the  choice  of  (6.12)  for  6V 
the  principle  of  virtual  work  does  indeed  imply  the  reduced  equation  of 
equilibrium,  and  so  (6.20)  contains  the  associated  natural  boundary 
conditions. 


6.2  Classical  Boundary  Conditions  alon 


•leral  Curves 


As  it  stands,  expression  (6.20)  is  somewhat  unsatisfactory  in  that, 
while  the  reduced  equation  of  equilibrium  (5.42)  and  the  variation  in 

Ct  ft 

strain  energy  (6.12)  depend  only  on  the  symmetric  part  of  M  it  appears 

that  (6.20)  involves  Ma^  entirely.  We  now  show  that  this  dependence  is 

spurious  by  first  writing,  without  loss  of  generality,  the  antisymmetric 
ot8 

part  of  M  as 


[aB]  _  aB 


(6.22) 


It  then  follows  that 


(f  [(M^  n).  „  *  6r  -  n  *  6aJ  v  ds 

.  I  p  _  -  _c  a 


=  $  (M  n  •  6r) ,g  e6“  ds  , 


(6.23) 


which  by  (A. 9)  becomes 


(j>  n)jg  •  6r  -  n  •  Sag]  ds 

C 

=  -  j  (M  n  •  6r)  ,g  x6  ds 


(6.24) 


Since  x  is  the  unit  tangent  to  C  and  satisfies 


ds  ds  _a 


(6.25) 


we  have  as  a  consequence  of  M  n  •  6r  being  continuous  that 


|  [(MM  n)|  g  •  5r  -  n  •  <5ag]  ds 


~  (M  n  •  6r)  ds  =  0 


(6.26) 


and  hence  that  (6.20)  depends  only  on  the  symmetric  part  of  M 


OB 
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Consequently,  by  virtue  of  (5.40) ^  we  can  write 


(Q  •  6r  +  m  •  68)  ds  =  (I  [(M6“  g  +  Qa  +  C°)  •  6r  -  M8a  •  Sag!  va  ds  . 


(6.27) 


Because  m  has  no  normal  component,  see  (4.21)2,  we  can  assume 
without  loss  of  generality  that  the  contact  couple  m  also  lacks  a 
normal  component 


m  •  n  =  0 


(6.28) 


Hence,  we  car.  define  the  equivalent  vector  M: 


M  =  m  x  n 


m  =  n  *  m 


(6.29) 


from  which  it  follows  that 


fi  fi 

1  •  50  =  -  H  •  a  n  •  6a „  =  -  M  n  •  6a 

-  -  -  -  -P  ^  -P 


(6.30) 


Substituting  this  expression  in  (6.27),  we  obtain 


{[(MBa,g  +  Q“  ♦  Ca)va  -  Q]  •  6r  -  [M6a  va  -  M8]  n  •  6aB}  ds  =  0  . 


(6.31) 


The  variation  in  aB  in  the  above  equation  cannot  be  independent  of 
6r  because  aa  =  r,a  .  To  isolate  this  dependence,  we  first  observe  that 
by  the  Kirchhoff  relation  (6.3) ^  6n  cannot  have  a  normal  component  and 
hence  can  be  resolved  tangential  to  S  along  and  normal  to  the  curve  C: 


on  =  t  x  •  6n  +  v  v  •  6n 


(6.32) 
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The  conyonent  v  •  6n  represents  a  variational  rotation  of  the  normal  n 
about  the  tangent  vector  r  and  is  independent  of  Sr  in  the  sense  that 
v  •  6n  can  vary  arbitrarily  when  Sr  =  0.  Hence  we  write  this  component 


as  a  rotation: 


6<j>  =  v  •  5n 


The  conponent  t  •  6n  however  is  entirely  determined  by  Sr,  since 


,  r  „dr  dSr 

T*Sn  =  -n*S'  =  -n*  S-f1-  =  -  n  •  -z — 

ds  ds 


With  the  last  two  expressions  5n  becomes 


6n  =  6$  v  -  n  *  t  , 


and  so  by  the  Kirchhoff  relation  (6.3} ^  we  have 

£ 


”  *  6aB  =  T8  "  *  df  "  VB  6<? 


(6.33) 


(6-34) 


(6.35) 


(6.36) 


Substituting  this  result  into  the  boundary  integral  (6.31)  and  inte¬ 
grating  by  parts  along  C,  we  obtain 


([(Me“B  *  *  C»)  V.  -  Q  *  i  (MB“  vo  tb  -  M6  n  t6))  .  6r 


(6.37) 


+  [(M8a  va  -  M8)  vg3  S4)  ds  -  l  [(M6a  vq  -  M6  n)  t  ]  •  Sr  1  =  0, 

1  ~  '  s+ 

i 

wherein  the  last  sum  is  over  points  i  of  discontinuity  in  the  tangent  t, 
and  hence  in  the  normal  v,  and  possibly  in  the  applied  contact  couple 
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component  MS.  Using  the  arbitrariness  of  5r,  we  can  separate  the  last 
equation;  first  setting  6r  =  0,  we  get 


[0Ta  va  -  MP)  vg]  64,  ds  =  0 


(6.38) 


which  when  substituted  back  in  (6.37)  gives  us  for  arbitrary  6r 


f  f!CM6a,6  *  ?°  -  f)  -  Q]  .  ^  liifa  v0  -  M°  n)  t6]>  •  «r  ds 


(6.39) 


.  I  [,W(Me“  Mws)  -  {’(«*'  v'-’-  n  •  «r  -  0, 


wherein  we  use  superscripts  (+)  and  (-)  to  denote  the  positive  and 

* 

negative  one-sided  limits. 

Using  the  boundary  conditions  (6.38)  and  (6.39)  we  obtain  the 
classic  boundary  conditions  along  a  general  curve  C: 

1.  Clamped  or  cantileverred  edge,  characterized  by  the  position  r  and 
the  normal  n  being  fixed  along  the  boundary;  hence  along  any  portion 
of  the  boundary  that  is  clamped  (6.38)  and  (6.39)  are  automatically 
satisfied,  since  6r  =  0  and  6<fr  £  0,  and  we  have  along  a  clamped 


r  =  const. 


n  -  const. 


(6.40) 


lim  j$(s)  and  li®  6(s)>  c^-  p.79]. 
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2.  Hinged  or  simply  supported  edge ,  characterized  by  the  contact  couple 
M  vanishing  and  the  position  r  being  fixed  along  the  boundary;  hence 
along  any  hinged  portion  of  the  boundary  (6.39)  is  satisfied,  since 
6r  3  0.  Assuming  the  integrand  continuous,  (6.38)  will  be  satisfied 
for  arbitrary  641  only  if  the  integrand  vanishes,  so  that  along  a 
hinged  edge 


r  =  const. 


=  0 


(6.41) 


3.  Free  edge ,  characterized  by  the  contact  force  Q  and  contact  couple 
M  vanishing  along  the  boundary,  edge  displacements  and  rotations 
being  arbitrary.  Hence,  assuming  the  integrands  continuous,  (6.38) 
and  (6.39)  will  be  satisfied  along  smooth  portions  of  the  boundary 
only  if 


(M1 
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ca) 


+  ^  (j* 


V  =  °» 


M 


aB 


-  0 


(6.42) 


and  at  comers  where  t  is  discontinuous  only  if 


(vW  -  ,<->  ,W)  -  0  . 


(6.43) 


6.3  Boundary  Conditions  along  Coordinate  Curves 

Along  coordinate  curves,  where  either  S1  or  £2  is  a  constant,  the 
free  edge  conditions  and  to  a  smaller  extent  the  hinged  edge  conditions 
take  on  sinpler  forms.  Let  us  first  consider  a  =  const,  coordinate 
curve  as  a  boundary.  Along  such  a  curve,  the  curve  tangent  x  is  colinear 
with  a2  and  the  curve  normal  v  is  colinear  with  a1.  Since  x  and  v  are 
unit  vectors,  we  have  that 

x1  =  0  ,  x2  =  ±l//a22  ,  vj  =  ±l//a^  ,  v2  =  0  ,  (6.44) 
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with  the  upper  sign  for  the  case  when  v  and  a1  have  the  same  sense  and 
the  lower  when  opposite.  Moreover, 


ti  =  a12 


t2  =  a22  t< 


=  .2  L 


(6.45) 


Substituting  these  results  in  (6.41)  and  (6.42),  we  obtain  for  edge 
conditions  along  =  const: 

2' .  Hinged  edge  (51  =  const) 


r  =  const.  ,  M*11  =  0  , 


(6.46) 


3'  .  Free  edge  (S1  =  const) 


M*11,!  +  2  M*12,2  ♦  N*1  =  0  ,  M*11  =  0 


(6.47) 


wherein  the  starred  variables  defined  by  (5.44)  and  (5.46)  are  most 
conveniently  used.  It  is  interesting  to  observe  after  writing  the 
reduced  equation  of  motion  (5.45)  in  extended  form 


Y*o  v  +  y*i  u>  +  [a  •  (y*i  v  +  y*2  “)  n] , 


=  (M* 1 1 , j  +  2  M*12,2  +  N*1) ,!  *  (M*22,2  +  N*2),2  +  F*  , 


(6.48) 


that  precisely  those  variables  differentiated  with  respect  to  51  in  the 
reduced  equation  are  the  ones  that  vanish  along  the  free  edge. 
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Next  taking  a  g2  =  const,  coordinate  curve  as  a  boundary,  we  have 
t1  =  +l//an  ,  x2  =  0  ,  vx  =  0  ,  v2  =  ±1//  22  »  (6.49) 

a 


with  the  upper  sign  when  v  and  a2  have  the  same  sense  and  the  lower 
sign  otherwise.  Moreover, 


A  ^ 

ti  =  au  x1  ,  x2  =  a21  x1  ,  -=-  =  x1 - .  (6.50) 

65  H1 

Substituting  these  results  in  (6.41)  and  (6.42),  we  obtain 
2".  Hinged  edge  (g2  =  const) 

r  =  const.  ,  M*22  =  0  .  (6.51) 

3".  Free  edge  (g2  =  constj 

M*22,2  +  2  M*12,!  +  N*2  =  0  ,  M*22  =  0  .  (6.52) 


Rearranging  the  reduced  equation  of  motion  (6.48): 


Y*o  v  +  y*i  w  +  [a  •  (y*i  v  +  y*2  J>)  a],o 


=  (M*22,2  +  2  M*12,!  +  N*2),2  +  (M*11,!  +  N*1),!  +  F*  , 


(5.53) 


we  again  observe  that  the  variables  differentiated  with  respect  to  g2 
vanish  along  the  free  edge. 
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const. 


When  a  boundary  coincides  with  the  intersection  of  a  = 
and  a  C2  =  const. coordinate  curve,  then  the  free  edge  condition  at 
comers  (6.43)  holds,  giving  us  in  addition  to  conditions  (6.47)  and 
(6.52)  the  condition 

3'".  Free  comer  (C1  =  const  §  C2  =  const.) 


M*12  =  0 


(6.54) 


6.4  Symmetry  Boundary  Conditions 

A  symmetry  edge  or  boundary  is  not  a  physical  boundary  of  the  shell, 
but  is  a  middle  surface  curve  that  always  lies  in  a  plane  about  which 
the  shell  deforms  symmetrically.  Hence,  initially  and  for  all.  subsequent 
positions,  this  plane  will  be  a  symmetry  plane  of  the  shell.  It  is  con¬ 
venient  to  take  the  symmetry  boundary  as  a  coordinate  curve  and  to  have 
the  remaining  coordinate  curves  symmetrically  distributed  with  respect 
to  the  symmetry  plane.  Also  the  external  force  F  and  external  couple  c 
must  be  symmetrically  distributed. 

Let  us  first  take  the  symmetry  edge  to  be  the  C1  =  C1  coordinate 

o 

curve.  Also  let  the  symmetry  plane  be  located  a  perpendicular  distance 
Kj  from  the  origin,  and  oriented  relative  to  the  constant  orthonormal 
basis  k^  (i  =  1,2,3)  such  that  kj  is  normal  and  hence  k^  (a  =  2,3) 
parallel  to  the  plane  as  shown  in  Figure  6.1.  It  then  follows  that 
when  C1  =  4  ,  the  position  vector  lies  on  the  symmetry  plane,  so  that 

f  (5*,  C2 j  t)  •  kj  =  Kj  (6.55) 

Moreover,  using  the  simplifying  notation 


6+  =  6  +  AC1,  C2), 


6_  -  &  (5q  -  AC 


1  >2 
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Figure  6.1  Deflection  pattern  of  the  middle  surface  of  a  shell 
deforming  symmetrically  about  a  symmetry  plane.  The  intersection  of 
the  synsnetry  plane  and  the  middle  surface  is  the  symmetry  boundary, 
for  which  material  coordinate  is  a  constant  51  =  5*.  The  synmetry 
plane  is  a  perpendicular  distance  K2  from  the  origin  0  and  normal 
to  the  vector  ki  of  ti.r  constant  orthonormal  basis  k-. 


--■-■■■.  . .. 


g 
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we  have,  as  a  consequence  of  the  coordinate  curves  being  symmetrically 
located  with  respect  to  the  symmetry  plane,  that 


(r+  +  r_)  •  ki  =  2  r  •  k;  =  2  Kj  ,  (r+  -  r_)  •  kQ  =  0  (a  =  2,3). 

(6.56) 

With  these  relations,  we  obtain  in  Appendix  B  the  symmetry  properties 
of  the  basis  vectors,  normal  vector,  metric,  etc.  and  combine  these 
results  with  the  velocity  relations,  strain  relations  and  resultant 
relations  co  obtain  the  symmetry  properties  of  the  variables  appearing 
in  the  reluced  equation  of  motion  (5.45): 


(?i 

-  aJ  j  • 

ki  =  0  , 

(a*  ♦  a*) 

•  ^  =  o  , 

+  a2)  • 

*1  =  o  , 

(a2  -  a2) 

•  ka  =  0  , 

(6.57) 

+  n  )  • 

ki  =  0  , 

(n+  -  n_) 

•  ka  =  0  , 

+  v_)  • 

ki  =  0  , 

-  y_)  •  ko  =  ° 

K 

+  0)  )  • 

k,  =  0  , 

(“♦ 

■  “.)  *  ka  =  0 

(6.58) 

(f; 

+  F*)  • 

ki  =  0  , 

(f: 

-  Fy)  •  ka  =  °  , 

M*11  =  M*11  , 

M*12  = 
+ 

-  M*12  , 

M*22  =  M*22 

q;11  »  q:11  , 

Q*12  = 

-  q:12  , 

q:22  =  q:2?  , 

(6.59) 

N*1  =  - 

n:1  , 

N^2  = 

N*2 
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Clearly,  the  above  symmetry  relations  hold  at  positions  symmetrically 
located  with  respect  to  the  boundary  and  are  useful  in  the  finite  differ¬ 
ence  formulation  of  the  reduced  equations  of  motion  along  symmetry 
boundaries.  Th^se  relations  can  be  also  be  used  to  obtain  the  conditions 
on  these  variables  and  their  derivatives  at  the  boundary,  for  observe 
that  as  AC1  -*•  0  ,  we  have  that  j$+  -*■  &  and  &  fa  .  When  we  use  this 

fact  and  note  that  aj  and  a1  at  the  boundary  are  colinear  with  kj  we 
find  that 

v1  =  0  ,  ui1  =  0  ,  F*1  =  0  ,  M*12  =  0  ,  Q*12  =  0  ,  N*1  =  0  .  (6.60) 

Next,  taking  the  symmetry  edge  to  be  the  C2  =  C2  coordinate  curve, 
we  locate  the  symmetry  plane  at  perpendicular  distance  K2  from  the 
origin  and  orient  it  relative  to  the  basis  k^  such  that  k2  is  normal  and 
hence  k^  (t  =  1,3)  parallel  to  the  plane.  With  the  notation 

<+  =  6  U1,  Z2q  +  AC2)  ,  =  i  (51,  5 1  -  AC2)  ,  {  =  t  (C1,  C2)  , 

we  then  have  that  the  position  vector  satisfies 

(r+  +  r_)  •  k2  =  2  r  *  k2  =  2  K2  ,  (rt  -  r_)  •  kT  =  0  (t  =  1,3).  (6.61) 

As  shown  in  Appendix  B,  with  these  relations  we  obtain  the  symmetry  of 
the  variables  appearing  in  (5.42): 


(a2  -  a2)  •  k2  =  0 


(a2  +  a2)  •  k  =  0 

-T 


(a*  +  a*)  •  k2  =  G 
(n+  +  n  )  •  k  =  0 


(a*  -  a*)  •  kT  =  0  ,  (6.62) 

(n  -  n  )  •  k.  =  0 
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"Tt*'*  >n*», iwmw 


(V+  ♦ 

v.)  *  k2 

=  0  , 

(v+  - 

v_)  •  kT  =  0 

y 

(oj+  + 

a)  )  •  k2 

=  o  , 

K  ~ 

t>  )  •  k  =0 
--  -T 

,  (6.63) 

(F*  + 

F*}  •  k2 

=  0  , 

?:  - 

F*)  •  kT  =  G 

y 

M*11 

=  M*U 

,  M*12  =  - 

M*12  , 

M*22  =  m*22 

y 

+ 

“ 

7  + 

— 

T  — 

„ 

q*h 

=  Q*11 

N*1 

+ 

,  q:12  =  - 

=  N*1 

q:12  , 

N*2 

+ 

q*22  —  q*22 

=  -  N*2 

>  (6.64) 

When  we  take  the  limit 

&2  and  a2  are  colineai 

;  as  AC2  -*•  o  , 

■  with  k2  ,  we 

using  the  fact  that  at 

find 

the  boundary 

v2  =  0  ,  ur  =  0  , 

F*2  =  0  ,  M*12  =  0  , 

q*12  =  0>  n*2 

=  0  .  (6.65) 

S9 


ga«a-a4i^,.v i»- 


7.  MATHEMATICAL  FORMULATION  OF  THE  REPSIL  CODE 


In  this  chapter  we  begin  by  specialising  the  general  shell  theory 
developed  in  the  previous  chapters  to  the  theory  on  which  the  finite 
difference  formulation  used  in  REPSIL  is  based.  This  specialization  is 
accomplished  in  two  steps:  first  the  thin  shell  approximation  is  applied 
to  the  general  equation  and  then  further  restrictions,  both  physical  and 
mathematical,  are  imposed  to  yield  the  special  theory  on  which  REPSIL  is 
based.  The  chapter  concludes  with  a  description  of  how  this  theory  is 
implemented  in  REPSIL  and  of  the  computational  procedure  used  to  obtain 
a  solution. 

7.1  Equations  of  Thin  Shell  Theory 

The  theory  developed  in  Chapters  2  through  6  assumes  only  that  the 

deformation  undergone  by  the  shell  satisfies  the  Kirchhoff  hypothesis. 

As  noted  on  page  29,  this  assumption  entails  a  mild  restriction:  the 

deformation  must  be  such  that  the  minimum  principle  radius  of  curvature 

is  greater  than  half  the  thickness;  i.e.  5-  <  R  .  .  Otherwise,  the 
b  2  mm 

position  of  particles  would  not  be  unique. 

We  now  impose  a  more  stringent  restriction  on  the  deformation^by 
requiring  R^.^  he  so  much  larger  than  h  that  terms  of  order  |s  b^|  and 
higher  can  be  neglected  and  dropped  in  comparison  to  unity.  In  other 
words,  we  assume  that  for  sufficiently  thin  shells  thest  higher  order 
terms  have  little  affect  on  the  deformation.  The  resulting  equations 
we  call  the  ii:in  shell  equations  and  the  resulting  theory  -  thin  shell 
theory.  The  theory  developed  here  differs  from  the  classical  theory  of 
thin  shells,  in  that  the  latter  neglects  terms  of  order  U  b^j  in  all  but 
the  expressions  for  the  strain  .  This  difference  will  be  pointed  out 
in  the  development. 


A  clear  and  straightforward  exposition  of  classical  thin  shell  theory 
using  the  Kirchh;ff  hypothesis  is  given  in  [13;  Ch.  2],  Direct 
comparison  of  tha  classical  equations  as  set  forth  in  this  reference 
and  the  equations  of  the  present  work  is  however  made  difficult  by  the 
unnecessary  use  of  principal  curvature  coordinates  in  the  former. 

PRECEDING  PAGE  BLANK 
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We  will  be  mainly  concerned  with  obtaining  the  thin  shell  versions 
of  the  equations  associated  with  the  reduced  equation  of  motion  (5.45), 
which  form  the  basis  of  the  REPSIL  formulation.  The  specialization  to 
thin  shells  of  equations  associated  with  other  forms  of  the  equations  of 
motion  is  however  straightforward. 

Relations  involving  the  geometry  of  the  middle  surface  will  remain 
unchanged.  Only  geometric  relations  off  the  middle  surface  will  be 
altered  by  the  thin  shell  assumption,  due  to  their  dependence  on  i ;  bg 
through  the  variable  Ug  being  made  linear.  Hence,  while  the  metric  a^g 
is  unchanged,  the  metric  g  R  given  by  (2.46)  becomes 


saB 


aB 


-  2  ?  b 


aB 


(7.1) 


It  follows  from  (3.12)  that  the  incremental  strain-displacement  relation 
(3.20)  also  becomes  linear  in 


AEaB  =  2  [(!a  *  A^B  +  !B  *  A“'a  '  A%  ’  A^B3 

"  2?(?  *  A“>aB  ‘  A?  *  !o,B  '  A?  ‘  A-’aB)]  * 


(7.2) 


For  computational  convenience,  tne  REPSIL  formulation,  rather  than  use 
the  last  form  of  the  strain-displacement  relation,  uses  the  equivalent 
but  shorter  form 


AE 


aB 


7  Uaa  •  Au,g  ♦  a8  •  Au,a 


A“’a  *  t">B) 


-  2?(n“ 


4%e  *  ^  • 


(7.3) 


The  surface  normal-displacement  relation  (3.30)  is  unaffected  by  the 
thin  shell  assumption. 


Since  all  the  resultant  relations  of  Chapter  4  are  obtained  by 
integration  through  the  thickness,  they  are  all  altered.  In  particular 
for  the  zeroth,  first  and  second  resui*  .  moments  of  mass  density 
(4.15),  we  now  have 
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■2 

^  =  Jh 


P  Cl  -  C  bY)(5)a  dc  (a  =  0,  1,  2)  , 


(7.4) 


which  follows  from  the  linearization  of  the  determinant  (2.47) Using 
these  relations,  the  momentum-velocity  relations  (4.14)  are  automatically 
appropriate  for  thin  shells.  The  resultant  relations  for  the  external 
force  and  couple  *4.36)  become 

h 

Fa  =  |2  [b“  -  25  bC“  bB)]  dc  +  [t“  +  t“  -  h  (b(“  tB)  -  bC“  tB))]  , 


2 

F  =  |  [1  -  5  bB]  b3dC  ♦  [t5+  +  t3_  -  ^-bB  (t f  -  1 1))  ,  (7.5) 

~2 
h 

c°  -  f  [b°  -  2c  b'»  b»)  «  dc  .  |  [t“  -  t“  -  h  (b<“  t«.  h<“  t«)]  . 


The  stress  resultant  and  couple  resultant  relations  (4.28)  assume  the  thin 
shell  forms: 


q“3  =  |  [0aB  -  25  b(B  o-)c]  d5  , 
Jh  ' 


Q°  =  f  U  -  c  bY)  05a  d5  , 

'h  Y 


(7.6j 


w“3  -  f  ;-“B  -  2C  b(S  5  d{  . 
Jh  7 
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However,  of  more  direct  concern  to  the  RErSIL  formulation  are  the  thin 
shell  forms  of  the  resultant  relations  for  the  modified  membrane  and 
bending  components: 

h 

Q°6  =  f 2  [1  -  ^  bY]  oaB  d?  , 
jh  Y 

"2 

(7,7) 

h 

5“6  =  f2  [Cl  -  t  b*)  c-6  -  5  b<°  „*»]  5  as  , 

jh  Y  Y 


which  follow  from  (5.43).  With  the  above  thin  shell  expressions  for 
Yo,  Yi,  Y2»  F®,  F,  Ca,  Q*6  and  MaB,  the  reduced  equation  of  motion  (5.41) 
or  (5.45)  is  appropriate  for  thin  shell  theory  without  need  of  change. 

As  mentioned  before,  classical  thin  shell  theory  goes  further  by 
dropping  the  linear  terms  in  (7.1),  (7.4)  -  (7.7),  giving  us  the 
approximations : 


gaB  aaB 


•  =  L  p  u 


)  dx,  (a  =  0,  1,  2)  , 


h 

~1 


h 

a 


f «  f/  M  ”7  T 

.a  t  ,  a  .  a  a  _  ,3,  3  3 

"  b  d?  +  t  +  t_  ,  t  =  b  d?  +  t  +  t  , 

J  h  Jh 


Q*6  .  QaB 


aS 


h 

~1 


(7.8) 


(2 


3a 


a  dc  ,  Q  =  dr  , 

A 

2 


C°  =  [  b“  ?  d?  +  |  (t®  -  t®) 


aB  «  r.aB  f  aB  „  , 


s  M 


j  o  ;  ac  . 
;h 
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Notice  that  the  distinction  between  components  Qa^  §  Ma^  and  the 

-aR  Act8 

modified  components  Q  6  M  disappears  completely  in  classical  thin 
shell  theory.  Moreover  notice  that  if  the  linear  terms  in  the  strain 
expressions  (7.2)  5  (7.3)  were  also  dropped,  then  the  strain  would  not 

take  into  account  bending  and  sc  would  be  appropriate  for  membrane 

* 

theory  .  In  any  event,  the  theory  forming  the  basis  for  the  REPSIL 
formulation  retains  for  the  sake  of  consistancy  terms  linear  in  £  b„ 

P 

in  the  membrane  and  bending  components. 


7.2  Restricted  Thin  Shell  Theory  Used  by  REPSIL 

The  equations  on  which  the  finite  difference  scheme  used  in  REPSIL 
is  based  are  simplified  forms  of  the  thin  shell  equations  just  developed. 
The  simplifications  are  of  two  types:  the  first  arise  from  the  specific 
type  of  problems  that  REPSIL  is  designed  to  treat  and  the  effects 
predominating  therein,  with  no  mathematical  approximation  being  used. 

The  second  follows  from  the  mathematical  need  to  make  the  equations 
simpler  to  solve,  and  involve  approximations.  Also,  we  adopt  an 
elastoplastic  constitutive  relation  between  stress  and  strain  to  make 
the  formulation  complete. 

REPSIL  is  formulated  to  treat  the  large  transient  response  of  thin 
shells  subjected  to  blast  loads.  For  this  reason  the  surface  tractions 
are  assumed  to  be  predominantly  due  to  blast  pressures,  with  the  effects 
of  viscous  drag  being  nil.  Moreover,  the  effect  of  gravicy  through  the 
body  forces  on  the  deformation  is  negligible  compared  to  the  blast 
pressure  effect.  Hence,  the  traction  and  body  force  fields  have  the 
special  forms: 


t+  =  ~  p+  "  ,  t_  =  P_  n  ,  b  =  0  , 


(7.9) 


Membrane  theory  is  treated  in  [12;  Sect.  14.3]  and  [13;  Ch.  4). 
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which  when  applied  to  the  components  of  F  and  C  as  given  by  (7.5) 
yield 

F  =  -  [P+  -  P_  -  |  (P+  +  PJ]  n  ,  C  =  0  .  (7.10) 

The  shell  is  assumed  to  be  initially  homogeneous,  so  that  the  initial 
mass  density  pq  is  constant.  Hence,  expressions  (7.4)  for  the  resultant 
moments  of  mass  density  can  be  integrated  for  their  initial  values, 
which  gives 

ro  =  P0h  .  r,  .  po  h3  ,  r2  =  ijpoh3.  C7.il) 


Consequently,  from  (5.20)  we  have 


Y*o  =  P0  h  Ah  ,  Y*i  =  -  J2  Po  h-1  A2  ,  Y*2  =  PQ  h3  A*  (7.12) 


As  mentioned  earlier,  the  first  resultant  moment  of  mass  density  is  a 
measure  of  the  distance  from  the  middle  surface  to  the  cer.cer  of  mass 
for  the  section.  It  is  clear  from  (7. 11)2  that  for  an  initially  uniform 
thin  shell,  the  center  of  mass  will  coincide  with  the  middle  surface 


when  the  initial  geometry  of  the  shell  is  such  that  the  mean  curvature 

B1  vanishes,  which  is  tantamount  to  the  principal  radii  of  curvature 
1 

being  equal  but  oppositely  directed.  It  then  follows  that  the  first 
resultant  moment  of  mass  density  will  vanish  for  all  subsequent  time, 
even  though  the  mean  curvature  b^  may  no  longer  vanish.  In  particular 
this  holds  true  for  initially  flat  shells. 

Expressions  (7.10)  -  (7.12)  are  exact  within  thin  shell  theory. 

For  the  sake  of  simplicity,  we  now  go  further  and  assume  that  the  first 
order  terms  can  also  be  neglected  in  these  expressions: 

F  =  -  [P  -P]n=-APn  ,  C  =  0  , 


r0  -  £>0  h  ,  Pi  =  0  .  f2  =  yy  pQ 


ft  1 


Y*0  =  P  h  A2  ,  Y*i  -  0  ,  Y*2  =  T7  h  » 

U  X  ~  o 
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giving  the  classical  thin  shell  approximation.  In  addition,  an 
approximation  is  made  that  gives  an  explicit  finite  difference  scheme. 
The  [7*2  aa  *  <I>  n],Q  term  in  the  reduced  equation  of  motion  (5.45)  is 
assumed  negligible  in  comparison  to  the  y*o  v  term.  Since  by  virtue  of 
(2.31)  a)  is  like  a  gradient  of  the  velocity  v,  it  follows  from  (7.13)_  _ 
that  this  approximation  is  equivalent  to  neglecting  |h  7  v|  in 
comparison  to  |vj.  In  other  words,  the  wave  lengths  of  acceleration 
waves  are  assumed  long  compared  to  the  shell  thickness.  With  the  last 
approximation  it  is  to  be  expected  that  the  solution  using  this  theory 
will  lose  accuracy  whenever  higher  harmonics  of  the  shell  are 
significantly  excited. 

The  above  specializations  of  the  thin  shell  equations  affect  only 
the  equations  of  motion.  In  particular,  the  reduced  equation  of  motion 
(5.45)  becomes 


with 


Y*0  v  =  .  +  N*“,  +  F* 

-  ~  a8  ~  a 


Y*0  =  pQ  h  ,  F*  =  -  AP  a'2  n 


(7.14) 


-aB 


„aS 


(7.15) 


M*~M  =  M*“w  n  ,  N*a  =  Q*a6aB  +  rg  a  M*By  n  . 


The  membrane  and  bending  components  are  not  specialized,  but  follow 
from  the  thin  shell  expressions  (7.7): 


=  a** 


j  (1  -  ?  bY)  oaS  d? 


(7.16) 


M*a6  =  a-  f  [(1-5  bY)  o°6  -  c  b(“  cyS>]  5  d5  . 


'_h 

■7 


Defining  the  zeroth,  first  and  second  moments  of  stress  as 
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caB  =  I 
a  Jh 


oaB  (Oa  d?  (a  =  0,  1,  2)  , 


the  last  two  relations  can  be  exprc-^ed  as: 

Q*“6  =  a4  [t*6  -  E“S]  , 

0  Y  1 

M*aB  =  a^  [ZaB  -  bY  ZaB  -  bCa  ZyB)] 
1  Y  2  Y  2 


(7.17) 


(7.18) 


In  order  to  complete  the  formulation  it  is  necessary  to  specify  in 
what  way  the  stress  is  related  to  the  strain.  Before  doing  this,  we 
first  characterize  the  stress  field  by  assuming  that  a  state  of  plane 
stress  exists  in  each  lamella  parallel  to  the  middle  surface  and  that 
the  shell  is  composed  of  isotropic  material.  However,  for  an  isotropic 
material  the  strain  components  given  by  the  Kirchhoff  hypothesis,  see 
equation  (3.6),  are  inconsistant  with  the  plane  stress  assumption,  since 
the  vanishing  normal  strain  component  will  usually  result  in  a  non-zero 
normal  stress  component.  This  situation  is  resolved  in  the  usual  way  by 
using  the  Kirchhoff  hypothesis  to  determine  only  the  tangential 
components  of  strain,  with  the  normal  components  following  from  the 
plane  stress  conditions. 


Presently,  the  REPSIL  formulation  employs  an  elastoplastic  stress- 
strain  law  that  is  linear  in  the  elastic  range  and  is  capable  of 

modelling  strain  hardening  and  strain  rate  effects  through  the  use  of 

* 

the  mechanical  sublayer  model  .  Using  the  above  assumptions  we  derive 
in  Appendix  C  an  expression  relating  the  incremental  change  in  the  stress 
components  to  the  components  of  the  strain  increment  and  the  present 
values  of  the  stress  components: 


ActT  = 


E 

1+v 


K 


v 

1-v 


AEyJ  - 


AX 


r  a 
°8  - 


I-^v 

3'il-vT 


aY] 

YJ 


(7.19) 


See  [2;  Sect.  5.4.2]  for  a  cor.cise  description  of  the  mechanical 
sublayer  concept. 
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wherein  the  mixed  components  are  with  respect  to  the  basis  g^: 


»\  - 


*  e“Y  ^ 


O  *1 

°8  =  %y  ° 


(7.20) 


The  factor  AX  is  non-negative  and  its  value  is  determined  from  the 
von  Mises-Hencky  yield  condition: 

,,  a,  _  a  8  1  ,  y>2  2  2  -  -  , 

v  8'  6  a  3  1  y'  3o  *  ' 


(7.21) 


with  oq  the  uniaxial  yield  stress.  When  AX  =  0  in  (7.19)  results  in  a 
state  of  stress  o^  +  AOg  outside  the  yield  surface  <{>(0“  +  Ao“)  >  0,  then 
the  proper  value  of  AX  is  the  smallest  positive  value  that  satisfies 
<j)(o“  +  Ao“)  =  0.  Otherwise,  AX  =  0  is  the  proper  value  since  it  gives  a 
state  of  stress  on  or  within  the  yield  surface  $>(Oq  +  Ao^)  <  0.  The 

P  P 

value  of  AX  is  thus  determined  by  (7.21)  and  the  stress  increments  (7.19) 
are  completely  specified.  For  a  strain  hardening  material  the  stress 
components  and  yield  stress  pertain  to  each  mechanical  sublayer.  A 
procedure  for  dealing  with  the  possibility  of  complex  AX  has  been 
developed  in  [4]  and  is  presently  incorporated  in  the  REPSIL  code. 


7.3  Description  of  the  REPSIL  Formulation 


We  give  a  description  of  the  computational  procedure  employed  in 

the  REPSIL  code,  including  an  outline  of  the  algorithm  used  to  advance 

★ 

the  solution  one  time  step  .  The  computational  procedure  is  similar  to 
that  used  in  PETROS  1  and  2,  so  that  the  following  description  will 
facilitate  the  comparison  of  these  codes  in  the  next  chapter. 

Basically  all  these  codes  solve  finite  difference  forms  of  partial 
differential  equations  that  characterize  the  transient  deformation  of  a 
thin  Kirchhcff  shell.  Finite  differencing  with  respect  to  the  material 
coordinates  results  in  the  dependent  variables  of  the  theory  being 
defined  on  a  two  dimensional  Lagrangian  mesh  or  grid  of  points.  At 


A  detailed  description  of  the  REPSIL  code  and  its  computational 
formulation  will  be  soon  forthcoming  in  a  user's  manual. 
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interior  mesh  points  and  at  those  along  symmetry  boundaries  central 
finite  difference  operators  are  employed,  while  along  fixed  boundaries 
(either  hinged  as  clamped)  forward  or  backward  difference  operators  are 
used,  with  all  operators  being  of  order  jAC^j^  .  At  present  the 
Lagrangian  mesh  has  a  rectangular  boundary,  which  means  that  the  REPSiL 
formulation  is  restricted  to  shells  having  simply  connected  middle 
surfaces  bounded  by  four  smooth  arcs. 

The  shell  is  conceptually  divided  by  a  number  of  lamellae 
parallel  to  the  middle  surface  into  K  layers  within  which  the  strain 
components  are  assumed  constant  at  their  mid- layer  values,  making  the 
stress  components  similarly  uniform.  The  reason  this  is  done,  despite 
the  fact  that  the  strain  increments  are  simple  linear  functions  ot  £ 
as  can  be  seen  from  (7.3),  is  that  plastic  flow  soon  makes  impossible 
the  simple  analytic  characterization  of  the  variation  of  the  stress 
components  through  the  thickness.  Thus,  the  components  of  stress  and 
strain  are  defined  at  stations  k  =  1,...K  through  the  thickness,  as  well 
as  at  each  point  of  the  Lagrangian  mesh.  This  dependence  on  k  will  not 
be  explicitly  indicated  in  what  follows. 

Replacement  of  time  derivatives  in  the  shell  equation  by  their 
equivalent  finite  difference  operators  results  in  an  explicit  finite 
difference  scheme  characterizing  the  ma7iner  ir.  which  che  dependent 
variables  of  the  shell  change  during  a  time  mt<- r^ai  At.  Hence  knowing 
the  values  of  the  shell  variable  at  time  t,  we  are  able  to  calculate 
their  values  at  the  time  t  +  At.  The  finite  time  derivatives  employed 
are  centered  and  of  order  | At | ^.  The  interval  At  does  not  change  with 
time,  but  is  fixed  at  a  value  for  which  the  explicit  finite  difference 
scheme  is  stable. 

Before  a  calculation  can  begin,  the  material  properties  and  the 
initial  geometry  of  the  shell  must  be  specified.  An  initial  velocity 
and/or  initial  pressure  distribution  starts  the  calculation,  giving  the 


A  finite  difference  formulation  of  the  free  edge  conditions  has  yet 
to  be  added  to  the  REPSIL  code. 
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initial  change  that  the  shell  variables  undergo.  After  that  a  sequence 
of  calculations  are  repeated  cyclically,  using  when  appropriate 
pressure-time  data,  to  generate  later  and  later  values  of  the  c lell 
variables  progressively  until  some  presciibed  maximum  time  is  reached, 
at  which  calculations  terminate. 

With  the  above  brief  description  of  the  finite  difference 

formulation  and  computational  procedure  employed  in  the  REPSIL  code  in 

mind,  we  are  ready  to  describe  the  algorithm  that  forms  the  heart  of  the 

* 

sequential  computations  advancing  the  solution  one  time  step  .  To  keep 

this  description  simple  we  shall  not  treat  a  strain  hardening  or  strain 

rate  sensitive  material,  although  as  mentioned  earlier  REPSIL  can  treat 

materials  exhibiting  such  phenomena.  The  algorithm  is  a  recipe  for 

solving  the  set  of  algebraic  equations  that  result  from  finite 

differencing  the  shell  equations.  For  our  purposes,  it  is  sufficient 

to  exhibit  the  finite  time  derivative  explicitly,  continuing  to 

symbolize  the  difference  operators  with  respect  to  £a  bv  partial 
** 

derivatives  . 

Let  us  assume  that  calculations  have  progressed  n-1  time  steps  to 
the  time  t  =  (n-1)  At  and  that  the  position  vector  r  and  normal  vector 

—  j* 

n  are  known  at  each  mesh  point  and  the  stress  components  a  are  known 
at  each  mesh  point  and  station  through  the  thickness.  Also  the 
displacement  increment  Au  for  the  time  interval  t~  -  .  1)  At  to 

t  =  n  At  is  assumed  to  have  just  been  calculated  at  each  mesh  point. 
Using  the  notation 


A  less  detailed  description  is  presented  in  a  recent  report 
[5;  Sect.  II]. 

ft 

The  piocise  finite  difference  expressions  will  be  in  the  forthcoming 
user's  manual. 
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6  =  i  [ (n-l)At]  ,  4  =  i  [nAt]  ,  i  -  i  [(n+l)At]  , 


H  =  i  -  f  .  H+  =  t"  -  i  ,* 

our  objective  is  to  outline  the  sequence  of  calculations  by  which  the 

new  values  r,  n,  aot5,  Au+  of  the  shell  variables  are  generated  from  the 
*•  ^  ^  _  — 0.8 

previously  known  values  r  ,  n  ,  o  ,  Au.  This  sequence  basically 
involves  five  calculations,  performed  in  the  following  order: 

7,3.1  Current  Geometry  Calculation.  The  current  position  vector 
r  is  found  from  (3.11): 


r  =  r  +  Au  . 


(7.22) 


Finite  differencing  the  current  position  vector  r  gives  the  variables 
characterizing  the  differential  geometry  of  the  middle  surface  (see  the 
equations  of  section  2.2): 

a  =  r,  ,  a  =  a  •  a.  , 

-a  -.  a  ’  a  B  -a  * 

a  =  an  a22  -  (ai2)2  ,  n  =  ai  *  a2/a2  , 

a31  =  a22/a  ,  a12  =  -  ai2/a  ,  a22  =  an/a  ,  (7.23) 


bae  "  !  '  r>ae  • 


b6  =  aBY  b  , 
a  ya  * 


a  ag 
a  =  a  a. 


r  y  =  aY  •  r 

a6  l  * 


7.3.2  Strain  Increment  Calculation.  In  order  to  determine  the 
strain  increments,  it  is  first  necessary  to  find  the  incremental  change 
in  the  surface  normal  by  finite  differencing  the  displacement 
increments  Au  and  using  (3,30)  and  (3.31): 


An  =  -  n  •  Au,  ,  An  =  [a 

a  «.  -a  -  _■ 


wt*  ^ 

aa  +  - —  n]  aat5  An  .  (7.24) 

~  l+n»n  ~ 


Notice  that  this  notation  is  in  accord  with  that  employed  earlier  in 
Sections  3,2  and  3.3. 
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The  strain  increments  through  the  thickness  at  stations  C(k),  k  -  1,...K, 
are  calculated  from  the  thin  shell  relation  (7.3): 

AE  =  \  (a  •  Au,„  +  a  *  Au.  -  Au.  •  Au, 

ap  2  _a  _  6  .6  -  a  „  a  -  p 

(7.25) 

"  5  Of  *  a6  +  *  !*aB}  * 

Notice  that  the  formulae  for  An  and  A£  „  are  at  least  linear  in  Au, 

.  ap 

assuring  numerical  accuracy  by  circumventing  the  need  to  take  differences 
of  nearly  equal  terms. 

7.3.3  Elastoplastic  Stress  Calculation.  The  stress  calculation 
uses  mixed  components  of  stress  and  strain.  Hence,  it  is  necessary  to 
know  the  metric  gQg  and  inverse  ga®  at  stations  c(k)  through  the 
thickness.  These  follow  from  the  thin  shell  relation  (7.1)  and  the 
equations  of  section  2.3: 

2 

8ag  ~  aag  “  2?  b^g  ,  8  =  f>ll  822  “  (812)  » 

(7.26) 

811  =  822/8  ,  812  =  -  812/g  ,  g22  =  gll/g  • 

With  these  relations  the  mired  components  of  the  strain  increment  and 
the  stress  are  obtained: 


(7.27) 


(7.28) 


(7.29) 


are  determined  and  used  to  find  the  current  values  of  the  stress 
components  using  (7.19): 


with 


C 


(7.30) 


=  Gn  - 


a  -a  .  a 
a8  =  0  8  +  A06  * 


(7.31) 


The  value  of  AX  depends  on  through  the  yield  condition  (7.21).  When 
T  B 


o®  is  on  or  within  the  yield  surface 

P 


T  T.  .  T  _  _ 

o“  o6  -  I-  (crT)2  **  I-  a2 
Bo  o  y  o  0 


(7.32) 


then  AX  =  0  and  current  stress  components  are  o^;  otherwise,  AX  is  the 
smallest  positive  value  giving  a  current  stress  on  the  yield  surface 


a  8  1  ,  y.2  2  2 

0  0  -  —  (o')  =  — 

B  a  .>  v  y  o 


(7.35) 


In  either  event  a  value  for  AX  and  hence  the  values  of  o“  are  determined. 

P 

Finally  the  contravariant  components  of  the  stress  are  found 


aB  _  ay  8 
0  =  g  °Y 


(7.54) 


When  the  value  of  AX  is  complex,  a  procedure  developed  by  HUFFINGTON  [4] 
is  used  to  determine  the  stress.  As  pointed  out  in  Appendix  C,  the 
stress  calculation  is  not  exactly  centered  with  respect  to  time  due  to 

-o  Co 

the  use  of  o  .  to  calculate  c„. 

P  P 


7.3.4  Membrane  and  Bending  Resultants  Calculation.  The  zeroth, 
first  and  second  moment  of  stress  are  determined  using  (7.17) 

zaB  =  £  oa8  (0«  ^  (a  =  0,  1,  2)  (7.35) 

k 

wherein  the  summation  sign  over  k  is  used  to  emphasize  that  numerical 
integration  through  the  shell  thickness  is  involved.  Kith  the  above 
values  of  stress  m-'-ants,  the  membrane  and  bending  components  follow 
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from  the  thin  shell  relations  (7.18): 


5*“6  =  [r“e  -  b>  r“6)  , 

M*“6  =  a1*  [e“6  -  bY  e“b  -  b(“  e~s))  . 
1  Y  2  Y  2 


(7.36) 


7.3.5  Equation  of  Motion  Calculation.  Using  the  just  determined 
values  of  Q*a8  and  M*“8  and  the  a  priori  specified  values  of  the 
pressure  differential  AP  between  the  shell  faces,  the  vectors 


N*a  =  Q*“8  ag  +  M*3y  n  , 


M*“8  =  M*a8  n  ,  F*  =  -  4P  n 


(7.37) 


are  determined  (cf.  (7.15)),  The  finite  difference  form  of  the 
equation  of  motion  (7.14)  with  the  above  Values  gives  the  displacement 
increment  Au+  for  the  next  time  increment: 


*  r°.„  *  r1  • 


wherein  y*o  is  the  initially  determined  time  constant 


(7.38) 


Y*0  =  PQ  h  A  • 


(7.39) 


This  calculation  completes  the  cycle,  the  current  values  r,  n, 

o“8,  Au+  having  been  determined  from  the  previous  values  r  ,  n  ,  a  a8, 

Au.  Hence,  the  solution  has  been  advanced  a  time  step  from  t  =  (n-1)  At 

to  t  =  nAt.  Repeating  this  sequence,  the  computer  code  generates  later 

and  later  values  of  r,  n,  oa8  ,  Au  until  some  prescribed  maximum  time 

t  =  NAt  is  reached, 
max 

7.5.6  Energy  Balance  and  Strain  Calculations.  The  sequence  of 
calculations  just  outlined  is  vital  for  the  advancement  of  the  solution. 
Concurrent  with  this  sequence  there  are  performed  calculations  that, 
although  not  essential  to  the  solution,  do  provide  information  about  the 
current  state  of  the  shell.  More  specifically,  as  a  global  measure  of 


the  deformation,  the  work  done  by  the  pressure,  the  kinetic  energy  and 
the  elastic  strain  energy  are  calculated  for  each  time  cycle.  The 
difference  between  the  pressure  work  and  the  sum  of  the  kinetic  and  strain 
energy  indicates  the  amount  of  energy  degradation  caused  by  plastic  flow. 
Also,  as  a  local  measure  of  the  deformation,  the  physical  components  of 
strain  at  selected  points  on  the  bounding  surfaces  of  the  shell  are 
calculated  for  each  time  cycle.  These  components  are  exact  for  large 
strains.  The  equations  employed  in  the  REPSIL  code  for  the  strain 
calculations  are  derived  in  Appendix  D  of  this  report.  The  derivation 
of  the  equations  used  for  the  energy  calculations  is  reserved  for  a 
subsequence  report.  The  forthcoming  user’s  manual  will  describe  how 
both  the  energy  balance  and  strain  equations  are  implemented  in  the 
REPSIL  code. 
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8.  COMPARISON  OF  THE  REPSIL,  PETROS  2  AND  PETROS  1  CODES 


As  mentioned  in  the  introductory  chapter,  the  REPSIL  and  PETROS  2 
codes  are  developments  of  the  PETROS  1  code,  both  correcting  and 
generalizing  the  PETROS  1  formulation.  In  this  chapter  we  compare  the 
theoretical  and  computational  capabilities  of  these  codes.  First  the 
theoretical  formulations  used  in  these  codes  are  compared  and  it  is 
shown  how  the  theory  on  which  PEiROS  1  is  based  has  been  generalized 
and  made  less  restrictive.  Next  the  computational  improvements  and 
increased  capabilities  of  REPSIL  and  PETROS  2  over  PETROS  1  are 
qualitatively  compared.  The  results  of  these  comparisons  are 
summarized  in  tables  and  the  implications  of  the  differences  between 
these  codes  are  discussed  briefly. 


8.1  Comparison  of  Theoretical  Formulations 

The  REPSIL,  PETROS  1  and  PETROS  2  codes  are  based  on  theories  that 
are  quite  similar.  All  assume  the  Kirchhoff  hypothesis  to  hold  and  all 
restrict  themselves  to  thin  shells.  Often  the  formulations  of  these 
codes  are  based  on  the  same  equations.  Essentially  all  employ  the  same 
formulae  to  describe  the  differential  geometry  of  the  middle  surface. 
They  differ  principally  in  the  equations  on  which  1°  the  strain 
increment  calculation,  2°  the  elastoplastic  stress  calculation,  3°  the 
membrane  and  bending  resultant  calculation  and  4°  the  associated 
equation  of  motion  calculation  are  based  (cf.  Sections  7.3.2  -  7.3.5). 
We  will  discuss  these  differences  in  this  order,  combining  the  last 
two. 

8.1.1  Strain  increment  Formulation.  The  various  strain- 
displacement  relations  employed  in  these  codes  have  as  common  bases 
the  Kirchhoff  hypothesis  and  thin  shell  assumption.  The  precise  forms 
of  the  relations  can  be  summarized  as  follows: 


REPSIL  (cf.  (7.24)  S  (7.25)): 


AEa8  =  2  (!a  *  A“’6  +  !e  *  Lxi> a  "  A!f’a  *  A“’5' 


■  *  A^’a8  +  Ae  *  *>cu3  ’ 

a3  .  . 

a  AnQ  Ang 

An  =  -  n~  •  Au,  ,  An  = - , 

ct  .  .’a  *  ,  - 

1  +  n  •  n 


(8.1) 


PETROS  2  [2;  Eq.  6.18,  S.23  5  5.29)]: 

AEa6  =  T  C?a  *  A“>S  +  !b  *  A%  ~  A“>a  *  A“>83 

-  C(n  •  Au,ag  ,  An  •  r,a6  -  An  *  Au,^)  ,C8>2) 

-y  n 

An  =  -  n  •  Au,  ,  An  =  1  -  (1  -  a  P  An  AnD)  2  , 
a  .  »’u  *  v  a  8 

PETROS  1  [1;  Eq.  3.46  8  3.40]: 


AEaB  =  2  [(^a  *  A“»8  +  ?8  *  A%3  +  l'(?a  *  ™»3  +  ffi  *  A?’a3 

An  =  -  n  •  Au,  -  n,  •  Au  ,  An  =  0  , 
u  _  _’a  _*a  _  * 


C«.3) 


with  An^  fi  An  being  the  components  of  An  for  all  three  formulations: 


An  =  An  a  +  An  n  . 
a  ^  _ 


(3,1) 


From  (3.16)7  it  is  clear  that  PETROS  2  uses  an  expression  for  AE^g 
that  is  exactly  equivalent  to  that  used  by  REPSIL.  Also,  both  use  the 
same  expression  for  Ana.  The  substitution  of  (3.16)2  it.  (8„i.),  gives  a 
quadratic  equation  on  An  of  which  (8.2)g  is  the  solution  under  the  mild 
condition  n  •  n  >  0,  so  that  both  REPSIL  and  PETROS  2  use  equivalent 
expressions  for  An.  Hence,  REPSIL  and  PETROS  2  have  equivalent  strain 
increment  formulations. 
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In  order  to  compare  the  expression  for  AE^g  used  in  PETROS  1,  we 
need  to  commute  the  order  of  differentiating  terms  in  (8.1)^  using 
(3.24)  to  obtain  the 


4E«e  *  2  u?a  ‘  se>6  *  !e  ’  5“>«  ‘  4iS  ‘  4“>  i> 

*  «;«  ■  4?>e  *  ?b  '  4;-i  *  •  4“>6  *  ?-e  •  4% 

1 

-  AnS(j  *  Au>6  -  An,g  •  Au,a)j 


(8.5) 


When  this  expression  and  the  REPS I L  expressions  for  the  components  of 
An  are  linearized  by  neglecting  products  of  gradients  of  the  displacement 
increments,  there  result  the  equations 

AE«xB  =  I  [(?a  *  A“'g  +  U  * 

+  5(!a  *  A?’B  +  tB  *  A?’a  +  ?»a  *  A“’g  +  ?>B  *  A“»a)]  » 


An  =  -  n  •  Au,  , 
a  -  „*a  * 


An  =  0  . 


(8.6) 


These  equations  would  coincide  with  the  PETROS  1  expressions  for  AE  ., 
An^  6  An  if  the  term  n,a  were  negligible,  which  means  that  the  middle 
surface  would  have  to  remain  virtually  flat.  Hence,  PETROS  1  uses  an 
approximation  to  the  exact  REPSIL  strain  increment  formulation  which  is 
valid  only  when  displacement  gradients  are  small  and  the  shell  is 
virtually  flat  throughout  the  deformation.  Consequently,  for  the 
particular  shell  geometry  treated  by  PETROS  1,  a  cylindrical  panel,  it 
cannot  be  expected  to  give  good  results. 

8.1.2  Elastoplastic  Stress  Formulation.  Both  the  REPSIL  and 
PETROS  2  stress  formulations  account  for  the  variation  in  the  metric 
through  the  thickness,  while  the  PETROS  1  formulation  disregards  this 
variation.  Specifically,  these  codes  use  the  following  formulae  for 
the  metric  and  inverse: 

*  - 

Notice  that  this  expression  is  equivalent  to  the  REPSIL  expression  for 

AEaB  in  the  partial  differential  sense,  but  not  the  finite  difference 

sense. 


REPSIL  (cf.  (7.26)) 


gaB  =  aag  ~  2C  baP  * 


g  =  gll  g22  ~  (gl2) ‘ 


g11  =  g22/g  ;  ,  g12  =  -  gl 2/ g  »  g22  =  gu/g 


(8.7) 


PETROS  7  [2;  ^  5.106  5  5.107] 


gag  =  aag  '  2C  bag 


ag  aB  _  .  ag 

g  =  a  +  2x.  u 


(3.8) 


PET.IOS  1  [1;  Sect.  5.1.2] 


■’ag  ^ag  ^ag 


ag  ag  _  -ag 
g  =  A  =6 


(8.9) 


The  loir-  a'  ror  g  .  used  by  REPSIL  ind  PETROS  2  are  identical  and  are 
based  *  -  t  lu'  .iiin  shell  assumption.  On  the  other  hand  REPSIL  employs 


the  exi : 


;t'  .iiin  shell  assumption.  On  the  other  hand  REPSIL  employs 
averse  to  g^,  while  PETROS  2  uses  an  approximate  thin 


shell  inverse.  PETROS  1  implicitly  assumes  that  the  metric  and  inverse 
«  ~e  equal  to  the  initial  middle  surface  metric  and  inverse,  which  for 
the  particular  initial  geometry  treated  (the  cylinder  section)  equal 
the  Krc.  ecker  delta.  Hence,  PETROS  1  becomes  more  restrictive  than 
either  REPSIL  and  PETROS  2  by  assuming  the  classical  thin  shell 
approximation  and  neglecting  any  variation  through  the  thickness 
(cf.  (7.8)j).  Moreover,  the  distortion  of  the  middle  surface  is 
neglected.  Thus,  PETROS  1  uses  an  approximation  to  the  thin  shell 
expressions  used  by  REPSIL  and  PETROS  2  for  the  metric  g^g  and  inverse 
ga^,  which  is  valid  for  extremely  thin  shells  that  suffer  negligible 
middle  surface  distortions. 

REPSIL  and  PETROS  2  employ  mixed  components  in  their  expressions 
for  determining  the  stress,  while  the  PETROS  1  expressions  do  not 
distinguish  covariant  from  contravariant  components  due  to  the  special 
form  of  its  metric.  However,  for  the  sake  of  comparison  we  shall  also 
write  the  °ETROS  1  expression  in  mixed  components,  summarizing  the 
expressions  used  by  these  codes  for  determining  the  stress  as  follows: 
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REPSIL  (cf.  (7.28)  -  (7.31),  (7.33)) 
T 


oa  =  a 

°B  8  1+v 


•°  +  J~  (AE“  +  ~  aE^  «“)  , 
8  1+v  ^  8  1-v  Y  8 


Ca  -a 
°8  =  °  8  "  3[H°  Y  UB  ’ 


l-2v  -Y 

a  o , 


(8.10) 


a  a  A,  a 
a8  =  °8  '  AX  °8  * 

PETROS  2  {2;  Eq.  5.68  -  5.75] 

T 


A  AX2  -  2B  AX  +  C  =  0 


•s  -  *  &  K  *  &  $  • 


Ca  Ta  l-2v  ^Y  *Y 

a8  =  f8  ~  Td^T  y  6b  * 


(8.11) 


T  C 

a  c  o 

°8  =  °B  '  AX  °8  * 


A  AX2  -  2B  AX  +  C  =  0  , 


PETROS  1  [1;  Sect.  4.3] 


oa  =  o*“  +  (AE“  +  AEj  6“)  , 
8  8  1+v  v  8  1-v  Y  p 


Ca  _  Ta  l-2v  ^Y 
°B  =  °8  "  3 (1-v)"  y  6B  * 


(8.12) 


-  2B  AX  +  C  =  0 


T 

VV  4i  °s  > 

where  mixed  components  ere  given  by  (7.27)  5  (7.54)  and  the  coefficient 

A,  B,  C  are  defined  uniformly  as 

C  C_  .  C  C 

.  a  8  1  a  _p 

A  =  o_  o  —  o  aR  , 

B  a  3  a  p 


C  T„  ,  C  T_ 

a  8  1  a  8 

B  =  °8  °a  ’  3  °a  °8  ' 

r  TaT8  lTa>  2  2 
C  =  °B  °a  '  3  °a  °8  "  5  °o  > 
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(8.13) 


MiattittB 


™  11  . . . ' 


in  all  the  above  formulae.  The  first  difference  to  be  noted  in  these 

formulations  is  that  while  REPSIL  and  PETROS  2  determine  AX  from  the 

exact  quadratic  formulae  thus  assuring  a  stress  on  the  yield  surface, 

p 

PETROS  1  uses  a  linear  approximation  to  determine  AX  giving  a  stress  not 
necessarily  on  the  yield  surface.  BUFFINGTON  [4;  Sect.  V]  has  assessed 
the  possible  errors  resulting  from  the  use  of  this  linear  approximation, 
indicating  that  in  a  typical  problem  involving  the  large  transient 
deformation  of  a  cylindrical  panel  on  the  average  a  discrepancy  of  50% 
in  the  stresses  can  occur  which  initially  can  go  as  high  as  210%. 

Hence,  we  can  conclude  that  the  use  of  a  linear  approximation  for  AX 
in  this  type  of  problem  is  certainly  inappropriate. 

A  second  not  so  obvious  difference  involves  the  stress  components 

C 

used  to  determine  corrector  stress  components  o“,  PETROS  1  and 
T  B  T 

PETROS  2  using  instead  of  c  °  as  P.EPSIL  does.  The  use  of  o”  is 

P  P  D 

rationalized  in  [2;  Sect.  5.4.3]  by  the  observation  that  a  real  value 
for  AX  is  more  likely  to  result  by  its  use.  However,  this  ad  hoc 
reasoning  becomes  much  less  persuasive  when  it  is  realized  that  the 
plastic  flow  rule  (C.10)  is  violated  by  this  choice,  because  o“  lies 

P 

outside  rather  than  on  the  yield  surface.  When  we  add  to  this  the 

fact  that  the  procedure  using  a~“  outlined  in  [4;  Sect.  IV]  and  used  in 

REPSIL  has  proved  in  practice  completely  adequate,  giving  real  AX  for 

T 

all  reasonable  deformations,  the  choice  of  become  still  less 

P 

justified. 

Lastly,  we  point  out  that  all  three  codes  have  the  capability  of 
treating  strain  rate  sensitive  and  strain  hardening  materials,  and  do 
this  in  more  or  less  the  same  manner  using  a  mechanical  sublayer  model 
for  strain  hardening  and  an  exponential  strain  rate  dependent  yield  stress 
[2;  Sect.  5.4.2]. 

8.1.3  Equation  of  Motion  and  Associated  Resultants  Formulation. 

The  various  forms  of  the  equations  of  motions  used  by  these  codes  are 
limited  by  the  common  assumptions  that 


C  =  0  ,  Yl  -  0  ,  Y2  =  0  • 


(8.14) 


Hence,  comparison  of  these  equations  will  be  within  the  limits  imposed 
by  these  assumptions.  The  equations  of  motion  and  the  associated 
membrane  and  bending  resultant  relations  used  in  these  codes  can  be 
summarized  as  follows: 

REPSIL  (cf.  (7.14)  -  (7.16)) 


Y*0  v  =  (M*p“  n),Ba  *  N*a,a  +  F*  , 


Y*o  =  P, 


h  k2  ,  N*  =  Q*a8  ag  +  rgYtt  M*8y  n  , 


•H 


aa&  (1  -  ?  bY)  d?  , 


(8.15) 


-  f2  l‘a8U  -  b(“  o*»]  c  dc  , 

J  h  '  ' 


PETROS  2  [2;  Eq.  2.67,  2.59,  2.83,  2.87,  2.88] 

Y*0  v  =  (M*Ba  n)  ♦  N*“,  +  F*  , 

-  d  _  a  _  a  _ 


*  r«„“  m*BY  "  > 


’■,  =  ,0“  •  r-«*  ?6*r6v 


Q*“B  .  zH  f  (6*  -  c  bB)  Cl  -  ?  b*)  d™  d?  , 


M*“6  -  a'S  /h  '  !  ijf  -  l  b°)  oY“  ?.  dC  , 


(8.16) 


the  current  Christoffel  symbols  can  be  used  in  the  PETROS  1  formulation, 
the  right  hand  side  of  (8.17)^  makes  no  restriction  as  to  the  smallness 
of  the  membrane  deformation.  Finally,  comparison  of  (8.17)4  ^  with 
(7.8) j.  ^  relations  (8.17)^  $  with  (7.8)^  ^  shows  that  the  PETROS  1 
resultant  relations  assume  the  classical  thin  shell  approximate jn, 
neglecting  terms  of  order  | 5  bg|. 

8.1.4  Summary  of  Comparison  and  Conclusions.  The  difference 
between  the  formulations  used  in  these  codes  are  summarized  in 
Table  8.1.  From  it  we  are  able  to  make  a  number  of  useful  conclusions 
about  the  kind  of  problems  in  which  these  codes  can  be  expected  to 
give  good  results.  First,  we  emphasize  that  all  the  formulations 
neglect  inertia  terms  containing  yj  and  Y2>  so  that  as  pointed  out  in 
Section  7.2  the  accuracy  of  solutions  will  diminish  whenever  harmonics, 
of  the  same  or  shorter  wavelength  than  the  thickness  are  present. 

REPSIL  and  PETROS  2  are  virtually  based  on  equivalent  formulations, 
both  being  able  to  treat  the  finite  deformation  of  thin  shells  for 
which  the  Kirchhoff  hypothesis  is  a  good  approximation.  Only  in  the 
characterization  of  the  flow  rule  do  they  markedly  differ,  the  REPSIL 
formulation  being  more  in  line  with  classical  plasticity,  while  the 
PETROS  2  foimulation  approaches  the  classical  case  only  as  the  strain 
increments  become  smaller. 

PETROS  1  is  not  as  general  as  the  other  codes,  which  is  to  be 
expected  since  it  is  the  progenitor  of  the  other  two.  The  linearized 
strain-displacement  relation  restricts  it  to  small  displacement 
increments  for  which  rotations  are  small.  Also,  the  improper 
characterization  of  the  curvature  effects  in  this  relation  make  the 
formulation  appropriate  only  for  virtually  flat  shells.  Due  to  the 
formulation  of  the  elastoplastic  stress  calculation  and  to  a  smaller 
extent  the  equation  of  motion,  the  membrane  deformations  must  be  small. 
In  addition  the  flow  rule  characterization  and  the  linearized  AX 
calculation  require  for  accuracy  that  the  strain  increment  be  small. 
Hence,  only  for  problems  involving  the  deformation  of  thin  shells  that 
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remain  virtually  flat  and  undergo  little  membrane  deformation  and  small 
incremental  displacements  can  the  PETROS  1  code  be  expected  on 
theoretical  grounds  to  have  the  accuracy  of  REPSIL  or  PETROS  2. 

An  assessment  of  the  relative  merit  of  PETROS  1  and  PETROS  2  in 
predicting  the  known  results  of  a  physical  experiment  is  attempted  in 
[2;  Sect.  6.3.2  6  6.3.3] .  However,  due  to  the  apparent  lack  of  precise 
control  of  the  edge  conditions  during  the  experiment  and  the  use  of  a 
rather  small  number  of  mesh  points  in  the  codes,  the  superiority  of 
PETROS  2  over  PETROS  1  was  not  definitively  proven. 


8.2  Comparison  of  Maior  Computational  Differences 


Similarities  in  their  theoretical  formulations  make  these  codes 
broadly  similar  computationally.  All  solve  finite  difference 
representations  of  the  particular  partial  differential  equations 
forming  the  bases  of  each  code.  They  all  attain  solutions  by  repeating 
a  sequence  of  calculations  which  advance  the  values  of  variables  one 
time  step,  similar  to  the  sequence  described  in  Section  7.5.  On  the 
other  hand,  differences  in  the  theoretical  equations  forming  their 
bases,  as  described  in  Section  8.1,  have  resulted  in  differences  in  the 
way  these  codes  are  computationally  formulated.  Differences  also  arise 
from  dissimilar  finite  difference  representations  of  the  same  formulae. 
Moreover,  there  are  differences  uue  to  different  computational 
capabilities.  We  shall  not  go  into  all  the  computational  differences 
between  these  codes,  but  only  describe  the  more  salient  ones. 

The  various  forms  of  the  equation  of  motion  used  by  these  codes 
give  rise  to  major  differences  in  the  finite  difference  representations 
of  the  second  partials  of  the  bending  components  in  these  equations. 
From  (8.16),  and  (8.17^  it  is  clear  that  in  PETROS  1  and  PETROS  2  the 
bending  components  occur  in  essentially  the  same  combination,  viz. 

(M^°,  n),  ,  while  (8.15)  gives  the  combination  (M^a  n),  for  REPSIL. 

ts  —  ct  l  «.  3  Ci 

For  a  typical  diagonal  component,  say  M11,  PETROS  1  uses  at  mesh  point 
(i,j)  the  obvious  finite  difference  representation 


] 


^'iJiiiiutiiciiaihiiiiiM^iiaaJi^tii.MLiiba  EflWttoflt6MiiafMaaft  haitmuttaittMiii  . . . . 


■ 


[(Mn,i  n),!],  ...  =  - - -  [M11,.  0  ...  -  M11,.  ..)  n,.  .  .. 

»u  (1,3)  4(  1)2  1  (1+2,3)  (1,3)-'  -(1+1,3) 


'  (MU(i,j)  "  MU(i-2,i))  U(i-l,j)]  * 


(8.19) 


This  formula  suffers  from  the  disadvantage  that  the  domain  of  mesh  points 
for  the  values  of  M11  skips  the  values  at  the  points  neighboring  (i,j) 
and  extends  over  an  interval  spanning  five  points.  In  other  words,  the 
domain  of  M11  is  too  coarse  and  too  extended.  PETROS  2  circumvents  this 
situation  by  the  judicious  use  of  forward  and  backward  difference 


operators : 


[(Mn,i  n).l](iJ)  -  ~ f(3Mll(i+i,j)  -  4Mll(i,j)  +  ““(i-lj)5  -(i+l,j) 


(8.20) 


+  (3M11,.  .  -  4M11,.  -s  +  M11,.  ,  -0  n,.  .  ..] 

(i-1,3j  (1,3)  (i+l,3)  ~(i-l,3) 


resulting  in  a  refined  and  compact  domain,  in  which  only  the  values  at 
(i,j)  and  its  two  nearest  neighbors  are  required.  Due  to  the  different 
form  of  the  equation  of  motion  used  in  REPSIL,  this  problem  is  not 
even  encountered,  for  the  straight  forwarl  use  of  the  finite  difference 
operator  for  second  derivatives  gives 


[(M11  n),n](iJ)  =  [Mn(i+1>j)  -(i+l,j) 


(8.21) 


-  2M11,.  ...  n,.  ...  *  M11,.  .  ...  nf.  .  -J  . 
(1,3)  -(1,3)  U-1,3)  -(i-l, 3) 


Not  only  does  this  formulation  result  in  a  refined  and  compact  domain 
for  diagonal  bending  components,  but  also  a  refined  domain  for  n  by 
including  its  value  at  (i,j). 


The  various  definitions  of  the  membrane  and  bending  components 


lead  to  differences  in  the  storage  requirements  of  arrays.  From 
(8.17)4  5  it  is  clear  that  as  used  in  PETROS  1  and  M°S  are 
symmetric  2*2  matrices.  However,  as  shown  in  the  previous  section. 


the  effects  of  the  middle  surface  curvature  have  been  ignored  through 
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a 


the  neglect  of  s  bg  terms.  When  this  curvature  is  taken  into 

consideration,  the  direct  generalization  of  Q°^  and  M°^  results  in 

definitions  similar  to  those  of  PETROS  2,  as  given  by  (8.16).  Clearly 

a8  u8 

the  symmetry  property  is  lost;  Q*  and  M*  are  nonsymmetric  2  *  2 

matrices.  Hence  the  number  of  independent  elements  has  been  raised  from 

six  to  eight.  Since  the  equations  of  motion  necessitate  the  finite 

differencing  of  these  components,  their  values  are  most  conveniently 

stored  as  arrays.  This  means  that  in  generalizing  from  PETROS  1  to 

PETROS  2  it  has  become  necessary  to  store  two  additional  arrays.  On 

the  other  hand,  the  generalization  used  in  REPSIL  to  account  for 

*  a8  *  ct B 

curvature  effects  retains  a  symmetric  formulation  with  Q*  and  M* 
as  defined  by  (8.15)4  g,  both  symmetric  2x2  matrices.  Hence  the 
REPSIL  generalization  of  the  PETROS  1  formulation  does  not  necessitate 
the  storage  of  two  additional  arrays  for  the  membrane  and  bending 
components. 

Also  the  REPSIL  and  PETROS  2  formulations  use  the  starred 
variables  and  Q*a^,  M*a^  rather  than  the  variables  Qa^,  MaS 

as  does  PETROS  1.  This  results  in  the  reduced  use  of  the  Christoffel 
symbol  in  the  equations  of  motion,  since  partial  rather  thrri  covariant 
derivatives  predominate,  making  REPSIL  and  PETROS  2  computationally 
simpler  than  PETROS  1, 

As  already  mentioned,  these  codes  treat  similar  analytic 
relations  in  computationally  different  ways.  The  numerical  integration 
of  the  stress  components  and  their  moments  through  the  shell  thickness 
for  the  membrane  and  bending  components  is  accomplished  differently. 
REPSIL  takes  straightforward  sums  of  the  values  of  the  stress  and  its 
moment  at  the  middle-point  of  each  layer,  as  does  PETROS  1.  PETROS  2 
uses  Gaussian  quadrature  for  this  integration.  In  [2;  Sect.  6.2.2]  the 
superiority  of  Gaussian  quadrature  over  middle-point  integration  is 
demonstrated  in  a  small  deflection  elastic  problem  for  which  Gaussian 
quadrature  models  the  deformation  through  the  thickness  exactly. 

However,  it  does  not  necessarily  follow  that  this  superiority  would  be 
so  pronounced  in  a  large  deflection  problem,  wherein  gross  plastic 
deformation  is  present. 


. 


Also,  these  codes  use  different  finite  difference  operators  for 
the  same  boundary  conditions.  PETROS  1  approximates  the  two  boundary 
conditions  it  treats,  the  symmetry  edge  and  clamped  edge,  by  central 
difference  operators  of  order  jA?  |  .  However,  as  pointed  out  in  [3], 
the  clamped  edge  condition  as  formulated  in  the  code  itself  is  improper. 
Early  versions  of  REPSIL  employed  the  same  boundary  conditions  using 
the  same  operators,  but  with  the  clamped  edge  condition  correctly 
fo — ‘ulated.  The  present  version  of  REPSIL  retains  the  central  difference 
operator  for  the  symmetry  edge,  but  uses  forward  and  backward  difference 
oper -.tors  of  order  |A£  |  for  the  clamped  edge.  In  addition  REPSIL 
t reals  the  hinged  edge  using  forward  and  backward  operators  of  order 
|A£a  .  PETROS  2  treats  the  symmetry  edge,  clamped  edge  and  hinged 
edge  conditions  more  or  less  in  the  same  way  as  REPSIL.  It  also  has 
the  capability  of  dealing  with  the  free  edge  and  cyclic  edge  conditions, 
using  central  difference  operators  for  both  [2;  Sect.  4.3].  Except  for 
the  free  edge  difference  operator  which  is  of  order  |A£a|,  all  the 
operators  used  in  PETROS  2  are  of  order  jAf;a|^.  As  of  now  it  is  not 
clear  whether  central  or  forward  and  backward  operators  are  preferable 
and  more  work  needs  to  be  done  here. 

The  comparison  of  the  difference  operators  used  at  the  boundaries 
points  out  some  differences  in  the  computational  capabilities  of  these 
codes.  These  codes  also  have  different  capabilities  as  to 
accommodating  various  initial  shell  geometries.  As  already  noted, 

PETROS  1  can  only  treat  initially  cylindrical  panels.  REPSIL  and 
PFfROS  2  have  been  formulated  so  as  to  accept  any  initial  geometry, 
provided  the  middle  surface  is  simply  connected  and  bounded  by  four 
smooth  edges.  This  versatility  is  achieved  by  confining  the 
specification  of  the  initial  shape  of  the  middle  surface  to  one  simple 
subprogram,  thus  making  the  remainder  of  the  codes  independent  of  initia’ 
geometry.  At  present  REPSIL  has  initial  geometry  subprograms  for  a  flat 
plate,  cylinder  and  cylindrical  ogive.  PETROS  2  has  subprograms  for  a 
flat  plate,  cylinder  and  cone.  Also  through  the  manipulation  of  values 
at  the  boundaries  to  eliminate  the  dependence  on  one  of  shell  parameters 
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£  ,  PETROS  2  can  treat  certain  symmetric  deformations  of  similarly 
symmetric  shells;  namely,  axisymmetric  deformations  and  plane  strain 
deformations  [2;  Sect.  4.3.7  §  4.3.8]. 

Both  REPSIL  and  PETROS  2  are  capable  of  handling  more  general  load 
situations  than  PETROS  1.  PETROS  1  can  only  handle  an  initial  impulse 
velocity  loading.  REPSIL  and  PETROS  2  can  also  accept  a  time  and 
space  varying  pressure  loading.  This  loading  can  be  specified  either 
analytically  or  numerically. 

The  results  of  this  section  are  conveniently  summarized  in 
Table  8.1.  It  shows  that  computationally  REPSIL  and  PETROS  2  are 
more  or  less  equal  and  both  are  more  sophisticated  than  PETROS  1. 

REPSIL  has  a  slight  advantage  over  PETROS  2  in  the  use  of  symmetric 
membrane  and  bending  components,  while  PETROS  2  uses  the  superior 
Gaussian  quadrature  to  integrate  through  the  thickness.  As  far  as 
capabilities  are  concerned,  PETROS  2  has  the  most  options  available 
and  REPSIL  has  more  than  PETROS  1.  As  is  to  be  expected,  REPSIL  and 
PETROS  2  are  comparable  and  both  are  superior  to  their  older 
predecessor  PETROS  1. 
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APPENDIX  A.  SURFACE  INTEGRAL  THEOREMS 

We  derive  in  this  appendix  the  surface  integral  counterparts  of  two 
volume  integral  theorems  used  in  continuum  mechanics.  First  we  derive 
Reynold*s  transport  theorem,  which  characterizes  how  the  operations  of 
time  differentiation  and  space  integration  commute  for  a  time  dependent 
surface.  Let  S(t)  be  a  time  dependent  surface  in  3-space  generated  as 
the  image  of  the  region  O  in  2-space  under  the  mapping 


r  =  r  U“,t) 


E.a  t  n 


(A.  1) 


and  let  £  (£a,t)  be  a  time  dependent  function  defined  on  0.  The  integral 
of  over  S(t)  can  be  written  as 


jj  6  da  =  jj  ^  a*5  dC1  d£2 


(A. 2) 


where  a  is  the  determinant  (2.7).  By  virtue  of  (2.28),  the  time  deriva¬ 
tive  of  (A. 2)  becomes 


^  ^  da  =  |j  (j  +  7  •  v  ^)  a:  d?1  d?2 


(A.  3) 


See  BUCK  [16,  Eq.  7-3]. 


w* £Ssk£&*- 


■  '-^<»4^r-.^-»Si»^iA^^»  _s*>s 


Applying  (A. 2)  in  reverse  to  (A. 3),  we  obtain  Reynold's  transport 
theorem  for  surfaces: 


jj  <5  da  =  ||  (jj  +  V  •  v  <$)  da  .  (A. 4) 


Next  we  obtain  the  divergence  theorem,  which  converts  surface 

integrals  into  line  integrals  along  the  surface  boundaries ,  by  special- 

* 

izing  Stoke's  theorem  to  surface  vectors.  For  a  continuous  vector 
function  ^  defined  on  a  smooth  surface  S  bounded  by  the  piecewise  smooth 
curve  C,  Stoke’s  theorem  reads 


V*^‘nda  =  |iS*Tds  , 


(A.5) 


with  t  the  unit  tangent  to  the  curve  C.  The  positive  sense  of  t  along 
C  is  determined  as  usual  by  the  right-hand-screw  relative  to  ■'•he  surface 
normal  n  (see  Figure  5.1).  The  following  relations  will  hold  between 
n,  x  and  the  exterior  unit  normal  v  to  C: 


x  =  n  *  v  ,  v  =  x  x  n  ,  n  =  v  x  x  . 


(A.  6) 


Since  t  and  v  are  normal  to  n,  a  and  a  can  be  used  as  bases  to  resolve 
-  -  -a 

these  vectors 


a  a 

t  =  x  a  =  t  a 
_  _a  a  _ 


a  a 

v  =  v  a  =  v  a 
_a  a  _ 


(A.  7) 


For  a  proof  of  Stoke’s  theorem  see,  for  example,  KELLOGG  [1/;  p.  89]. 
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and  conversely  x,  v  can  be  used  as  a  basis  for  aQ  and  a  : 


a  a  a 

a  =  x  x  +  v  v  ,  a=TT+VT 
.a  a  _  a  .,  5 


Moreover,  the  components  of  x  and  v  satisfy  the  relations 


CA.8) 


a  ga 
x  =  e  v. 


v  =  e  0  x 
a  aB 


CA-9) 


If  we  let  ^  be  a  vector  tangent  to  the  surface  S,  so  that 


we  find  that 


V  x  j  .  n  =  n  •  a“  x  aP),a  =  e“e  ^ 


Stoke' s  theorem  then  becomes 


|j  e“B  <01  a  **  =  <f  <8  xB  ds 


Finally,  setting 


a  aB  / 
9  =  e  <P 


(A.  10) 


and  substituting  in  (A.  10),  we  obtain,  in  view  of  (A. 9),  the  divergence 
theorem  for  surfaces: 


jj  9°(  a  da  =  I  <?“  va  ds 


(A. 11) 


v 


APPENDIX  B.  SYMMETRY  BOUNDARY  FORMULATION 

A  symmetry  plane  is  a  plane  about  which  the  shell  is  always  sym¬ 
metrically  positioned.  Hence  initially  the  shell  is  symmetrically 
distributed  with  respect  to  the  symmetry  plane  and  subsequently  deforms 
in  a  symmetric  manner  about  this  plane.  Consequently,  the  body  force  b 
and  surface  tractions  t  are  required  to  be  symmetrically  distributed 
with  respect  to  the  symmetry  plane.  The  intersection  of  the  middle 
surface  with  the  symmetry  plane  gives  the  symmetry  edge  or  boundary. 

For  convenience,  we  take  the  symmetry  boundary  to  be  a  coordinate  curve 
and  require  the  remaining  coordinate  curves  to  be  symmetrically  dis¬ 
tributed  relative  to  the  symmetry  plane. 

We  first  let  C1  =  C1  be  the  symmetry  boundary.  The  symmetry 
plane  is  located  a  perpendicular  distance  Kj  from  the  origin  and  oriented 
relative  to  the  fixed  orthonormal  basis  (i  =  1,2,3)  such  that  k}  is 
normal  and  hence  k^  (a  =  2,3)  parallel,  as  shown  in  Figure  6.1. 

The  position  vector  to  any  point  on  the  symmetry  edge  by  definition 
satisfies 


r  (C1,  C2,  t)  •  kj 


(B.l) 


Also,  using  the  symplifying  notation 

£+  =  6  (C1  +  H1,  £2) ,  £_  =  £  U1  -  AC1,  C2),  £  =  £  U1,  C2)  , 

we  obtain  as  a  consequence  of  the  symmetric  distribution  of  coordinate 
curves  about  the  symmetry  plane 

[5+  *  r_]  •  k,  •  2  r  -  k,  =  2Kj,  [r+  -  rj  •  k,  =  0  (a  =  2,3).  (6.2) 
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The  symmetric  distribution  of  the  body  force  b  and  surface  tractions  t 
is  equivalent  to  the  conditions 


(i)+  '  b(D-  1 

b(o)+  b(o)- 

1 

t — s 
H 
v— ' 

•M 

1 

II 

+ 

r~\ 

H 

w 

t(o)+  "  t(a)- 

(B.4) 


The  symmetry  properties  of  all  other  variables  follow  from  (B.3)  and 
(B.4).  In  particular,  from  (4.33)  and  (4.34)  we  have 


F  =  -  F 

(D+  (D- 


F(o)+  F0»- 


C(1U  "  "  C(l)-  » 


C(o)+  C(o). 


(B.5) 


Before  continuing  with  determining  the  symmetry  properties  of  the 
remaining  variables,  we  must  first  indicate  how  the  symmetry  properties 
of  variables  change  under  differentiation.  The  above  relations  clearly 
show  that  the  conponents  of  vectors  with  respect  to  the  basis  are  of 
two  types:  symmetric  and  antisymmetric,  satisfying  respectively  the 
typical  relations 


~  — 
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It  follows  from  the  definition  of  the  derivative  that  differentiation 
with  respect  to  C1  will  change  the  symmetry  type  of  a  variable  while 
differentiation  with  respect  to  5  2  will  not,  so  that 


K  =  i. 


^*2+  ~  ”  ^’1-  *  <5 >2+  ^ *2-  * 


»  6 >2+  ”  6 >2_ 

With  these  properties,  we  obtain 

properties  of  the  basis  : 

on  differentiating  (B.3) 

al(l)+  =  al(l)- 

al(a)+  =  ‘  al(o)-  » 

a2(l)+  =  ’  a2(l)-  , 

a2(o)+  =  a2(o)- 

(B.6) 


from  which  the  symmetry  properties  of  the  metric  a^  follow 


a,,.  -  a,,  ,  a,,.  =  -  a,0  ,  a22+  -  a^ 


11+  11-  *  12+  12- 


(B.7) 


Using  (2.4)  and  (2.10)  we  find  the  symmetry  properties  of  n  and  a  : 


“(«♦  '  '  "(l)-  ’  n(a)+  =  "(o)-  * 


1  1 

a(D+  a(D- 


2  _  2 
a(D+  '  a(D- 


(o)+ 

2 


=  -  a 


1 


CO-  , 


(B  -  8) 


S(a)+  =  a(a)- 
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Differentiation  of  (B.6)  gives  us 


al,l(l)+  '  al,l(l)-  ’  al,l(o)+  3l,l(a)-  * 


al>2(l)+  al,2(l)- 


al  ,2 (o)  +  "  '  *1,2(0)-  , 


(B.9) 


a2,2(l)+  '  a2,2(l)-  ’  a2,2(a)+  =  a2,2(o)- 


Confcining  the  last  two  results ,  we  obtain  the  symmetry  properties  of 


b  .  and  r  Y  : 
aB  otB 


bll+  bll-  1  b12+  "  b12-  *  b22+  b22- 


rni  -  -  rni  >  ri2!  -  ri2?  ’  r22!  •  -  r22-  •  'B-10) 


r2-r2  r2  r2  r2-r2 
11+11-  ’  12+  ~  12-  ’  22+22- 


Tbe  definition  (3.14)  of  the  incremental  strain  AE  shows  that  it  has 

Otp 

the  same  symmetry  properties  as  a^  and  b^^ 

AE11+  =  AE11-  ’  AE12+  =  "  AC12-  ,  AE22+  =  AE22-  ’  C6-11) 


from  which  it  follows  that  the  stress  a  also  has  similar  symmetry 
properties: 


11  11 
a  =o 
+• 


22  22 
a  =a 
+ 


(B. 12) 


see  Appendix  C. 
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The  symmetry  properties  of  the  membrane  and  bending  components  we 
infer  from  their  definitions  (5,43),  from  which  it  follows  that 


Q1.1  »  Qn  ,  i\2  -  -  Q1.2  ,  ?}-$}  , 

if}  =  if}  ,  Hi2  =  -  M12  ,  Sf  =  if}  . 


(B.  13) 


Also  conbining  CB- 5) ^  4  with  (B.8)^  5  we  have 


c*  =  -  c1  ,  c2  =  C2 


(B.  14) 


Hence,  by  definitions  (5.40),  (5.44)  and  (5.46),  we  obtain 


Q*11 

x+ 

=  q:11  , 

Q*12 

=  -  i*12 

,  if2  =  Q*22  , 

A  11 
M* 

«  M*U  , 

m;12 

-  -  M*12 

,  M^2  =  M^22  , 

(B. 15) 

"1  ~  1 
N*  =  -  N* 


n*2  =  n*2 


The  symmetry  properties  of  the  velocities  v  and  u>  easily  follow  from 
taking  the  time  derivatives  of  (B.3)  and  (B.8). 


V(o)+  =  v(a). 


%)  +  "  to(o)- 


(B- 16) 


To  facilitate  obtaining  the  symmetry  conditions  when  the  symmetry 
boundary  is  the  £2  =  £2  coordinate  curve,  we  note  that  the  symmetry 
characteristics  of  all  the  above  variables  are  related  to  the  number 
of  times  the  index  1  appears:  a  variable  is  symmetric  when  the  index  1 
occurs  an  even  number  of  time  and  antisymmetric  when  it  occurs  an  odd 
number.  This  is  a  result  not  only  of  the  symmetry  edge  being  a  C1  =  £* 
coordinate  curve,  but  also  of  the  symmetry  plane  being  normal  to  basis 
vector  kj.  Hence,  for  the  case  where  5  2  =  5  2  is  the  symmetry  edge, 
the  symmetry  plane  is  positioned  a  perpendicular  distance  K2  from  the 
origin  and  oriented  relative  to  the  fixed  basis  k^  such  that  k2  is 
normal  and  hence  k^  (x  =  1,3)  parallel-  By  analogy  with  the  previous 
case,  we  now  have  that  variables  in  which  the  index  2  occurs  an  even 
number  of  times  are  symmetric  and  those  in  which  it  occurs  an  odd  number 
antisymmetric.  Thus,  we  have: 


r(2)+  +  r(2)- 


=  2  K„ 


r(T)+  “  r(-r)~  *  =  (B.17) 


F(2)+ 

‘  F(2)- 

F(t)+  : 

C(2)+  = 

"  C(2)~ 

CW  : 

al  (2)+  = 

"  aK2)-  ’ 

^(T)^ 

a2  (2)+  = 

a2  (2)  - 

a2  (t  )  + 

1 

1 

1 

a(2)+  = 

■  a(2) - 

^T)- 

2 

2 

2 

a(2)+ 

a(2)- 

a(T)  + 

"(2)+  = 

"  "(2)- 

n(T)+ 

(B. 18) 


(B* 19) 
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APPENDIX  C.  ELASTOPLASTIC  STRESS  FORMULATION 

We  present  here  a  brief  derivation  of  the  formulation  forming  the 
basis  of  the  stress  calculation  used  in  REPSIL.  The  derivation  is 
similar  to  that  given  by  HUFFINGTON  [4;  Sect. II],  but  generalized  to 
general  non-orthonormal  bases.  The  material  composing  the  shell  is 
assumed  to  be  isotropic.  Also,  the  surface  tractions  are  assumed  so 
small  compared  to  the  internal  stress  field  and  the  shell  assumed  so 
thin  that  the  shear  and  normal  stress  components  in  the  direction  of 
the  shell  normal  can  be  neglected.  Hence,  using  components  of  the 
stress  a  relative  to  the  basis  g^,  we  have 


? 


I. 


(C.l) 


| 

§ 

s- 

S’ 


This  gives  us  a  state  of  plane  stress  in  each  lamella  parallel  to  the 

middle  surface.  Using  mixed  components  of  stress,  the  nonzero  components 
a  * 

of  stress  are  Og. 

The  material  is  assumed  to  satisfy  a  linear  incremental  stress  - 
strain  relation  in  the  elastic  range.  In  the  plastic  range,  the  stress 
components  are  delimited  by  the  von  Mises-Hencky  yield  condition 
[18;  p.20],  which  for  the  case  of  plane  stress  in  mixed  component  fcrm 
is 


(C.2) 


with  a  the  uniaxial  yield  stress.  Thus,  the  stress  is  restricted 
o 


Raising  and  lowering  of  indicies  is  accomplished  through  the  use  of 

metric  g  0  and  not  a 
saB  aB 

This  treatment  applies  to  the  mechanical  sublayer  model  of  strain 
hardening  by  considering  the  stress  Og  and  yield  stress  a  to  be  those 
in  an  individual  sublayer.  Moreover,  strain  rate  sensitivity  is  also 
included  by  making  oq  a  function  of  the  rate  of  strain.  Cf.  [2;  Sect. 
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to  values  within  the  yield  surface  <  0)  or  on  the  yield  surface 
(<}>  =  0).  Plastic  straining  only  occurs  when  the  stress  has  values  on 
the  yield  surface. 

We  formulate  an  incremental  theory  of  plastic  behavior,  assuming 

that  the  incremental  strain  components  AE^*  are  the  sums  of  corresponding 

B  E 

components  of  the  elastic  strain  increment  AE0  and  the  plastic  strain 
P  6 

increment  At  : 

p 


E  P 

._a  .ra  ..-a 
^8  =  ^6  *  ^6 


(C.3) 


The  elastic  strain  increment  components  are  related  to  the  stress  incre¬ 
ment  components  through  an  isotropic  Hooke's  law: 


A_a  1+v  .a  v  .  y  ta 
4EB  ’  —  %  -  E  4%  66  • 


(C.4) 


Solving  (C.4)  for  the  stress  components  and  using  (C.3)  we  find 


40b  =  K  ‘  *  iVs  • 


(C.5) 


Hence,  if  strain  and  hence  stress  increments  occur  in  the  time  interval 
[t  -  At,  t],  we  have 


Og(t)  =  (t  -  At)  +  Aa“ 


(C.6) 


Clearly,  the  state  of  stress  at  the  time  t  -  At  can  be  presumed  to  be 
an  allowable  state,  so  t'  at 


*(o“(t  -  At))  <  0  . 


(C.7) 


**■•  .W.  1  V^^tS^-  -  T-^-  Z 


Tentatively  setting  to  zero  the  plastic  strain  increments  in  (C.5),  if 
the  state  of  stress  at  the  time  t  is  allowable,  so  that 


<Koe(0)  <  0  , 

p 


(C.  8) 


then  indeed  the  change  in  the  stress  state  has  occurred  elastically 
and  the  stress  at  the  time  t  is 


a«(t)=a“(t-At)+1Ar[AE“+I^6“  AE^] 


(C.9) 


If  however  the  stress  at  the  time  t  is  not  allowable: 

<KOg(t))  >  0  , 

then  the  plastic  strain  increments  cannot  be  assumed  to  vanish,  but 
must  be  determined  from  the  flow  rule.  Using  the  yield  condition  (C.2) 
as  a  plastic  potential  (see  [18;  p  51]),  we  obtain  for  the  flow  rule 


,,-a  3<f>  ,r  _  ,  a  1  ,a  y.  ~ 

6  ^  C"e  *  5  6  V 

a 


(C.10) 


When  we  write  the  flow  rule  in  a  form  more  appropriate  to  our  formu¬ 
lation  using  finite  increments: 


,ra  _  .  CX  1  ,a  Y.  ~ 

=  2  (°B  -  5  6  6  V  AX 


(C.ll) 


the  question  immediately  arises  as  to  what  values  of  the  stress  com¬ 
ponents  to  use  in  this  formula.  Since  the  strain  increments  are 
associated  with  the  time  internal  [t  -  At,  t] ,  ideally,  values  of  the 

components  of  stress  should  be  for  some  intermediate  time,  say 
a  l 

°g  (t  -  2  At) ,  in  order  that  the  stress  be  properly  centered.  However, 
such  intermediate  values  are  undetermined.  Moreover,  using  the  average 
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of  the  values  at  t  -  At  and  t,  vis. 


1 


j  [o“(t  -  At)  +  o“(t)]  ,  would  not 


only  lead  to  computational  difficulties  since  the  values  o0(t)  are  as 

p 


yet  unknown,  but  would  be  theoretically  wrong  since  the  average  would 
not  give  a  state  of  stress  on  the  yield  surface  due  to  the  strict  con¬ 
vexity  of  this  surface.  For  these  reasons  we  are  persuaded  to  use  the 
values  at  the  time  t  -  At  for  the  stress  components  in  (C.ll).  With 
this  understanding,  we  substitute (C. 11)  in  (C.5)  to  obtain 


A  a„  = 


E  [AEa  +  v  6u  aEy 
1+v  1  B  1-v  fi  y 


-  2  {ag(t  -  At) 


l-2v 


3(l-v)  B  Y 


6®  aY(t  -  At)}  AX]  ,  (C.  12) 


or  using  the  abbreviations 


.  a  C  r,r“  v  -a  .rli 

Aa„  =  - -  [Ac  +  - - o.  At  ] 

B  1+v  *•  B  1-v  B  Y 


C 


_a _ a  l-2v  xa  _y 

°B  B  ‘  3(J-v)  6B  Y 


(C. 13) 


2  £  - 
A  _  *  1 


AX  =  AX 
1+v 


obtain  the  compact  relation 


Aa“  =  Aag  -  AX  c“(t  -  At) 


(C. 14) 


We  now  can  proceed  to  determine  the  values  of  the  stress  components  at 
the  time  t  by  substituting 


MO  =  °«(t  -  At)  +  AOp  (AX) 


(C- 15) 
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^ i 


in  the  yield  condition  (C.2)  resulting  in  a  quadratic  expression  in  AX 
(notice  the  dependence  of  Act®  on  AX  in  (C.15),  which  we  solve  for  the 
smallest  positive  value  of  AX  satisfying 


<HOg(t))  =  0 


This  value  of  AX  will  satisfy  the  quadratic  equation 


where 


A  AX  -  2B  AX  +  C  =  0 


c  c  c 

A  =  o“(t  -  At)  o£(t  -  At)  -  ±  [o“(t  -  At))2  , 


(C. 16) 


(C- 17) 


T  C 

B  =  o“(t)  o*(t  -  At) 


C  =  <HOg(t)) 


5  “sct  '  «> 


•C-18) 


°g(t)  =  Og(t  -  At)  +  Ac“ 


Details  of  how  the  value  of  AX  is  determined  from  (C.17)  as  well  as  a 
procedure  for  dealing  with  the  case  of  complex  AX  can  be  found  in 
[4;  Sect.  IV). 


S 

y 


We  can  summarize  the  procedure  for  determining  the  values  of  the 

T 

stress  components  as  follows:  first  a  trial  stress  o  (t)  is  determined 
from  the  previous  stress  °g(t  -  At)  and  the  strain  increment  using 
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3 


»*■  ~  ,  - 


«5 


T 

(C.13),  and  (C.18)..  If  the  trial  stress  oD(t)  satisfies  the  yield 
14  p 

condition: 

T 

<J>(Og(t)<0  ,  (C.  19) 

then  the  trial  stress  is  the  stress  at  the  time  t: 

T 

o“(t)  =  a“(t)  .  (C.20) 


If  the  trial  stress  does  not  satisfy  the  yield  condition: 

T 

<Ka“(t))  >  0  ,  (C.21) 


then  solving  (C.17)  for  a  proper  value  of  AX,  the  stress  at  the  time  t 
is  given  by 


C 

a“(t)  =  a“(t)  -  AX  a“(t  -  At)  ,  (C.22) 


C 

wherein  a0(t  -  At)  is  determined  from  the  stress  at  the  time  t  -  At 

P 

using  CC- 1 3) 2 - 
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APPENDIX  D.  PHYSICAL  COMPONENTS  OF  SURFACE  STRAIN 

In  this  appendix  we  derive  the  expressions  used  in  the  REPSIL  code 
for  the  physical  components  of  strain  on  the  bounding  surfaces  of  a 
shell.  These  expressions  hold  for  a  thin  Kirchhoff  shell  and  simulate 
the  readings  of  strain  gages  bonded  to  the  shell.  They  are  exact  with 
no  restrictions  as  to  the  smallness  of  strain  components  being  imposed 
in  the  derivation. 

in  Chapter  3  we  derived  expressions  (3.2)  and  (3.3)  for  the  current 
and  initial  differential  arc  lengths  between  neighboring  shell  particles. 
For  particles  on  the  bounding  surfaces  of  the  shell,  c  =  and  hence 
dt  =  0,  so  that  these  expressions  become 


ds2  =  ga6  d5“  d56 


d52  =  GaB  d5“  d£* 


(D.l) 


The  strain  E  that  a  fibre  connecting  neighboring  particles  on  the 
surface  undergoes  is  by  definition 


E  =  —  -  1 
dS  1  ’ 


(D.2) 


This  gives  us 


C£  +  -  SaB  T“ 


(D.  3) 


,  a  _  d£ 

where  x  =  ^g—  are 

basis,  since  clearly 


where  x  5  -r|—  are  the  components  of  a  unit  vector  x  in  the  initial 


G  _  Ta  tB  =  1 
ap 


(D.4) 


Substituting  definition  (3.5)  in  (D.3)  and  taking  account  of  (D.4),  we 
relate  the  physical  strain  E  to  the  tensor  components  of  strain  E^: 


(E  +  l)2  =  2EqB  t°  t8  *  1  . 


(D.5) 


Since  t  is  a  unit  vector  in  the  plane  of  the  basis  Gq  ,  it  is 
completely  specified  by  the  angle  0^  it  makes  with  or,  equally  well, 
the  angle  0O  it  makes  with  G_.  Hence,  we  can  write  its  components  as 
follows 


x1  =  G1  •  T  =  Jip-  sin  9.  ,  T2  2  G2  •  T  =  v^22  sin  8,  .  (D.6) 


Clearly,  0^  and  9 are  not  independent,  but  are  related  through  the 


metric  G 

a0 


cos  (Gj  +  e2)  = 


/G11  G22 


,  sin  (0J  +  02)  !c  G  ,  CD. 7) 


'11  22 


with  G  the  determinant  of  the  metric: 


G  =  G11  G22  ~  'G12^ 


(D.8) 


Using  (D.7)2,  we  can  rewrite  (D.6): 


1  _  _L 


sin  9, 


2  1 


sin  9, 


/GjY  sin  (0X  +  02) 


vG22  sin  (0J  ♦  02) 


(D-9) 


* 


r 


i^ 


Substitution  of  the  last  expressions  in  (D . 5)  gives  us  the  strain  E  as 


a  function  of  the  angle  the  fibre  makes  with  the  basis  G^: 


(e  + 1  r  = 


2  ,  ? 
-  [z^  sin  $2 


sin  (6 x  +  02) 


(D. 10) 


.  2 


+  2y  sin  6^  sin  02  +  e2  sin  0^]  +  1  , 


where  we  use  the  abbreviations 


’31 


'1  G 


Y  = 


12 


11 


SgTTg 


Z2  G 


"22 


11  12 


22 


(D. 11) 


Successively  setting  0^  and  @2  equal  to  zero,  we  obtain  the  strain  in 


the  G^  and  G2  directions: 


1  =  A  +  2  eJl  -  1  ,  E2  =  A  +  2  e2  -  1  . 


(D.  12) 


Since  0^  ai  d  ©2  are  dependent,  we  may  eliminate  either  one  from 
L)  through  the  use  of  ( 
the  trigonometric  identity 


CD- 11)  through  the  use  of  (D.7).  We  proceed  to  eliminate  62  by  dividing 


sin  0?  =  sin  (0^  +  ©2)  cos  0,  -  cos  •''.9,  +  0,,)  sin  0 


1  2 


1 


CD.  1 3) 


by  sin  (9^  +  9?)  and  using  (D.7)  with  the  abbreviations 

G, 


o  = 


12 


A  = 


'AlS? 


h  -  62 


(D. 14) 
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55. 


*Zi,  .  «>.  *„,^*t-  ,iy.-s^Sv/<  05“  -w-Sr  rft/Y'S**1*'  »;•*  »'-*^ 


-  *,  *  -Jjt- ■ 


to  obtain 


sin  0. 


6(6^6) 


sin  (61  +  82) 


cos  0j  -  5  4  sin  6^ 


(D. 15) 


Using  (D.7)  and  D.14)  again,  we  also  have 


sin  9, 


“(6i>6)  5  sr„- (if  *  o,t  ~ 4  si"  ei 


(D. 16) 


On  substituting  the  last  two  expressions  in  (D.10)  and  solving  for  E, 


we  obtain  the  strain  for  a  fibre  in  the  direction: 


E  =  [2(S2  e2  +  2  a  B  Y  +  a2  e2)  +  1]^  -  1  . 


(D. 17) 


Equations  (D.12)  and  (D.17)  are  the  ones  used  in  REPSIL  to  compute  the 
physical  strains  on  the  bounding  surfaces  of  the  shell. 


Lastly,  we  show  that  when  the  basis  is  initially  orthogonal,  so 
that  Gj2  =  0,  and  the  strain  components  are  small,  e^,  e2  and  y  become 
the  usual  small  strain  components;  for  clearly  we  then  have  from  (D.12) 


E1  *  el  »  E2  *  e2 


(D. 18) 


iso  from  (D.l)  and  (D.ll)2  we  have 


- hi - =  1  (  hi-*22-  j  cos  (i  -  2r)  , 

2>?I7S7  l  G11  °22  I 


(D.  19) 


whe.  ;  2f  is  the  change  the  angle  between  the  basis  vectors  experience 


is  going  from  Gq  to  g^.  Using  (D. 11)  and  (D.12),  the  last  expression 


(D.20) 


i 

t 

! 

| 


i 

t 


becomes 

Y  =  \  (1  ♦  Ej)  (1  +  E2)  sin  2r  , 
which  for  small  strains  is  approximated  by 


Y  *  r  .  (D.21) 

Hence,  for  small  strains  and  approximate  the  strain  undergone  in 
the  G^  and  direction  while  y  approximate  half  the  angle  change  under¬ 
gone  between  the  directions  G^  and  G2;  these  are  the  usual  small  strain 
definitions  for  the  strain  components  in  the  plane. 


